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ABSTRACT 


This  thesis  is  concerned  with  approximation  by 
rational  functions . 

Basic  Chebyshev  concepts  regarding  polynomial 
approximation  are  reviewed.  Theorems  concerning  the 
approximation  of  a  continuous  function  by  a  rational 
function  are  presented  along  with  proofs.  Then  basic 
Chebyshev  concepts  are  employed  in  an  attempt  to  outline 
a  constructive  theory  of  approximation  by  rational 
functions . 


Some  algorithms  for  obtaining  rational 
interpolatory  functions  and  rational  approximations  have 


been  included. 
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CHAPTER  I 


INTRODUCTION 


Section  1.1  The  Nature  of  Approximation 

In  many  of  the  problems  which  arise  in  numerical 
analysis.,  the  mathematician  is  given  certain  information 
concerning  a  particular  function  that,  for  example,  may  be 
known  to  be  -  or  assumed  to  be  -  continuous.  He  is  then 
required  to  educe  additional  or  improved  information  in  a 
form  which  is  appropriate  for  interpretation  in  terms  of 
numbers . 

A  technique  which  is  frequently  used  in  such 
cases  can  be  described,  in  general  terms,  as  follows.  A 
convenient  set  of  n+1  coordinate  functions,  say  cpq(x), 
cp^(x) ,  .  .  .  . ,  cpn(x),  is  first  selected.  Then  a  procedure 

is  invented  which  has  properties  such  that  it  will  yield 
the  desired  additional  information  easily,  and  (except 
for  inaccuracies  in  calculation)  exactly  if  f(x),  the 
function  to  be  approximated,  is  a  member  of  the  set, 

Sn-n  maybe  infinite,  of  the  functions  which  are  express¬ 
ible  exactly  as  linear  combinations  of  the  coordinate 
functions.  Here  we  are  assuming  that  an  appropriate 
criterion  for  approach,  such  as  the  Chebyshev  criterion 
(see  below  and  Chapter  2),  least  squares  criterion,  etc.. 
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has  been  chosen.  Different  criteria  will,  in  general, 
yield  different  approximations,  and  any  particular 
criterion  will  reproduce  some  properties  of  f(x)  more 
faithfully,  other  less  faithfully.  Some  loss  of  fidelity 
must  be  accepted  in  an  approximation;  but  -  as  mentioned 
above  -  the  selective  process  should  arrive  at  f(x)  if 
f(x)  is  in  S  . 

One  further  practical  point  is  recognized.  In 
general,  the  first  approximation  to  a  function  uses  a 
bare  minimum  of  information  and  thus  is  crude.  The  process 
must,  in  general  be  such  that  the  approximation  can  be 
improved  by  incorporating  additional  information.  It  is 
also  desirable  that  an  existing  approximation  should  be 
able  to  use  additional  information  concerning  f(x),  but 
not  used  by  the  selective  process,  for  estimating  the 
error  of  the  approximation. 

Obviously,  it  is  advantageous  to  choose  coordinate 
functions  whose  properties  are  simple  and  which  are  easily 
computed.  Since  polynomials  are  easily  evaluated  and  since 
their  integrals,  derivatives,  and  products  are  also  poly¬ 
nomials,  the  n+1  functions,  1,  x,  ...,  xn,  which  generate 
the  algebraic  polynomials  of  degree  ^  n ,  are  particularly 
appropriate  -  again,  n  may  be  infinite. 
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The  history  of  polynomial,  and  other  approximations 
spans  two  centuries  and  only  a  brief  outline  of  develop¬ 
ments  of  major  importance  will  be  given  here  -  not  neces¬ 
sarily  in  chronological  order. 

Section  1.2  The  Beginning  of  Approximation  and  Interpolation 

Interpolation  may  be  defined  as  any  method  for 
obtaining  numerical  values  of  a  function  which  uses  only 
certain  tabular  values  of  that  function. 

It  appears  as  though  Gregory  should  be  regarded 
as  the  father  of  interpolation  based  on  the  representation 
of  functions  by  power  polynomials.  Although  Newton’s 
interpolation  formulae  using  forward  and  backward  differences 
first  appeared  In  print  In  his  "Philosophiae  Naturales 
Principia  Mathematica" ,  Oxford,  1687  (Book  III,  Lemma  V), 
they  were  in  fact  discovered  by  Gregory  some  seventeen 
years  earlier.  Briggs  originated  the  use  of  tabular 
differences  of  higher  than  the  first  order  for  the  construc¬ 
tion  and  interpolation  of  tables  In  his  "Arlthmetica 
Logarithmia" ,  London,  1624,  and  later  in  "Trigonometric 
Britannica",  Gaudae,  1633)  which  he.  and  Gellibrand  wrote. 
Briggs  appears  to  have  preceded  Gregory  but  he  did  not 
develop  any  Interpolation  formulae  using  the  higher 
order  differences.  The  terms  "interpolatio"  and  "interpolare " 


4. 

in  their  present  sense  were  apparently  introduced  by 
Wallis  in  his  "Arithmetica  Inf initorum" ,  Oxford,  1659 . 

Newton,  in  his  "Principia"  (Book  III,  Lemma  V) 
gave  the  general  notation  of  divided  differences  while  the 
term  itself  appears  to  have  been  coined  by  DeMorgan  in 
his  "Differential  Calculus",  London,  1842.  In  addition, 
Newton  gave  us  two  interpolation  formulae,  commonly  known 
as  Sterling’s  and  Bessel's,  first  with  unequal,  then  with 
equal  argument  intervals,  in  "Methodus  Dif f erentialis " 
of  the  "Analysis  per  quantitatum.  Series,  Fluciones,  ac 
Diff erentias . . . London,  1711.  It  can  also  be  shown 
that  Gauss's  interpolation  formula  is  only  a  slight 
transformation  of  Sterling's  (Newton's)  formula  (see  p.66 
of  [18]). 

A  result  of  fundamental  importance  in  numerical 
work  is  the  Lagrange  interpolation  formula  which  was 
closely  anticipated  by  Euler  in  his  "Institutiones  Calculi 
Dif f erentialis " ,  Petrograd,  1755^  in  which  he  reached 
his  formulae  by  considering  series  of  which  a  certain 
order  of  difference  is  constant.  The  highest  polynomial, 
one  of  the  fifth  order,  that  Euler  gave  was  deduced  from 
Newton's  forward  difference  formula.  However,  the  Lagrange 
formulae  was  given  explicitly  by  Waring  in  "Phil.  Trans. 


: 
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Roy.  Soc."  69,  59  (1779)  and  was  subsequently  published 

by  Euler  in  his  "Opuscula  Analytica",  (Vol.  I),  Petrograd, 

1783^  in  the  form  appropriate  for  use  with  equidistant 

data.  Apparently  this  latter  publication  does  not  antedate 

Waring,  but  it  does  antedate  Lagrange,  who  did  not  publish 

the  formula  that  is  now  generally  associated  with  his 

✓ 

name  until  1795  in  "Jour,  de  l'Ecole  Polytechnique",  2, 

274.  The  logical  connection  between  Lagrange’s  (Waring* s) 
and  Newton’s  divided  difference  interpolation  formulae 
was  first  explicitly  pointed  out  by  Gauss. 

Another  central  difference  interpolation  formula, 
commonly  known  as  Everett's  formula,  was  published  in 
1901  by  Everett  in  "Jour.  Inst.  Act."  _35,  452.  However, 
the  formula  was  actually  first  derived  by  Laplace  in  his 
"Theorie  Analytique  des  Probability " ,  Paris.  This  formula, 
unlike  other  central-difference  formulae,  uses  only  even- 
order  central  differences,  and,  as  a  result,  if  only 
lower-order  differences  are  used,  requires  less  tabulation 
of  differences.  Another  formula,  usually  referred  to  as 
Steffenson’s  interpolation  formula,  was  credited  to 
Everett  by  Steffenson  who  called  it  Everett's  second 
formula  in  his  "Interpolation". 


It  can  be  shown  (see  p.l6  of  [ 18 ] )  that  the 


l 


forward-difference,  central-difference,  and  the  Lagrangian 
formulae,  when  the  abscissae  are  equally  spaced,  are 
really  different  aspects  of  the  same  process;  that  is, 
the  running  of  a  polynomial  of  degree  n-1  through  n  points. 
Hence  when  we  use  these  interpolation  formulae,  we  are 
assuming  that  the  interpolation  polynomial  and  the  function 
coincide  closely  enough  to  give  an  error  of  interpolation 
that  lies  within  prescribed  bounds.  However,  because 
the  higher-order  differences  involved  in  the  forward  and 
backward-difference  formulae  are  less  and  less  closely 
related  to  the  behaviour  of  the  function  near  the  desired 
point,  these  interpolation  formulae  are  not  as  accurate 
in  practice  as  other  difference  f ormalae .  Thus  the  forward 
and  backward-difference  formulae  should  only  be  used  in 
cases  where  it  is  impossible  to  use  other  interpolation 
formulae;  that  is,  at  the  ends  of  the  difference  table. 

Even  then,  it  is,  if  possible,  better  to  extend  the  differ¬ 
ence  table  and  use  a  central-difference  interpolation 
formula . 

Section  1.3  The  Modern  Theory  of  Approximation 

Approximation  maybe  described  as  the  act  of 
replacing  a  known,  untractable  function  by  a  tractable 
function  having  similar  characteristics. 
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Let  us,  at  this  point,  exemplify  how  the  criterion 
of  approximation  influences  the  nature  of  the  interpolatory 
approximation.  (All  the  approximations  considered  at  this 
point  are  polynomial  approximations) .  The  use  of  zeros 
of  the  Legendre  polynomials  in  approximation  minimizes  the 
root  mean  square  value  of  the  error  over  the  open  interval 
(-1,  l)  whereas  the  use  of  the  zeros  of  the  Chebyshev 
polynomials  minimizes  the  largest  of  the  absolute  values 
of  the  errors.  Furthermore,  the  errors  oscillate  uniformly 
over  (-1,  l),  while  in  the  first  case  they  tend  to  oscillate 
with  increasing  amplitude  towards  the  ends  of  the  interval. 
It  was  shown  by  Chebyshev  in  1853  that,  when  the  maximum- 
error  criterion  is  used,  the  choice  of  the  Chebyshev 
polynomials  is  the  best  possible  one.  In  the  series  of 
theorems  containing  the  above  result  he  gave  a  practical 
and  simple  construction  which  enables  one  to  find  excel¬ 
lent  approximations  under  mild  restrictions.  A  discussion 
of  Chebyshev  concepts  will  be  found  in  Chapter  II. 

A  natural  extension  of  polynomial  approximation 
is  approximation  by  rational  functions;  that  is,  by  the 
ratio  of  two  polynomials.  Cauchy,  in  his  "Cours  d’analys£ 
de  I’Ecole  polytechnique  ,  premier  Partie,  Note  V,  Paris, 
1821,  proposed  an  interpolation  formula  involving  the 
ratio  of  two  polynomials  of  different  orders.  Thus, 
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Cauchy  was  one  of  the  first  to  investigate  the  possibility 
of  interpolation  and  approximation  by  rational  functions. 

In  fact,  it  appears  to  be  very  difficult  to  find  any 
branch  of  analysis  upon  which  Cauchy  has  not  touched.  In 
1846,  E.  Brassinne  published  a  paper  in  "Journal  de 
Mathematiques  pures  et  appliquees",  Tome  XI,  Paris,  which 
indicated  the  variety  of  formulae  which  are  included  under 
Cauchy's  ratio  formula.  Theile,  in  his  "Interpolations- 
rechnung",  B.  G.  Teubner,  Leipzig,  1909*  published  his 
reciprocal-difference  algorithm.  This  algorithm  allows 
the  construction,  from  a  table  of  reciprocal  differences, 
of  a  rational  approximation  in  continued  fraction  form, 
which  is  the  easiest  and  most  efficient  form  for  evaluation 
of  rational  functions  of  interpolation  type.  The  algorithm 
will  be  given  in  Chapter  IV. 

Up  to  this  point,  we  have  considered  mainly 
interpolating  approximations.  This  type  of  approximation 
may  be  characterized  as  approximation  to  differentiable 
functions.  A  radically  different  method  of  approximation 
was  initiated  by  Weierstrass  when  he  proposed  to  consider 
approximations  to  continuous  functions  by  a  polynomial 
of  sufficiently  high  order.  Weierstrass  published  his 
theorem  on  polynomial  approximation  in  "Sitzungsber  Akad 
Berlin",  in  1885.  This  famous  theorem  states,  in  effect. 
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that  any  function  continuous  over  a  closed  interval  can 
he  uniformly  approximated  over  that  interval,  to  any 
prescribed  tolerance,  by  a  polynomial  of  sufficiently 
high  order.  Bernstein’s  theorem,  published  in  1912,  is 
a  significant  improvement  over  Weierstrass ' s  result, 
since  it  actually  constructs  a  sequence  of  polynomials 
such  that  the  approximation  improves  steadily  as  the  order 
of  the  approximating  polynomial  increases.  However, 
Bernstein’s  method  converges  too  slowly  to  be  of  much 
practical  value.  Today  there  is  much  interest  in  the 
Chebyshev  theorems  and  with  modern  high-speed  digital 
computers  we  are  increasingly  able  to  find  Chebyshev- 
approximants  to  continuous  functions. 

Section  1.4  Survey  of  the  Thesis 

The  remainder  of  this  thesis  will  be  concerned 
with  the  following  topics.  Chapter  II  consists  of  a 
presentation  of  the  Taylor  polynomials,  the  Lagrange 
formula,  the  Newton  polynomials  and  some  Chebyshev  concepts 
regarding  polynomial  approximations.  Chapter  III  contains 
an  extension  of  Chebyshev  concepts  to  rational  approxima¬ 
tions,  while  chapter  IV  contains  some  algorithms  for 
finding  rational  interpolants,  and  chapter  V  contains 
algorithms  that  give  -  in  most  cases  -  "best-fit"  (in  the 
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Chebyshev  sense)  rational  approximations.  Chapter  VI 
contains  a  summary  of  what  has  been  covered  in  this 
thesis  and  a  discussion  of  some  significant  unsolved 
problems  concerning  rational  approximations. 
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CHAPTER  II 

POLYNOMIAL  APPROXIMATION 


Section  2.1  Introduction 

This  chapter  contains  a  discussion  of  the  Taylor 
polynomials  and  the  development  and  a  discussion  of  the 
Lagrangian  formulae,  the  Newton  polynomials  and  some  basic 
Chebyshev  concepts  concerning  polynomial  approximation. 

The  Taylor  polynomials  are  included  because  they  are  among 
the  forms  most  widely  used  on  today’s  high  speed  computers. 

When  the  notation  f(x)  is  used  in  what  follows, 
it  is  meant  to  represent  an  arbitrary  continuous  function 
defined  in  a  finite  interval  a<x<b,  and  in  section  2.2, 
f(x)  is  differentiable  any  number  of  times  at  some  point  c 
contained  in  the  interval  of  Interest. 


Section  2.2  The  Taylor  Polynomial 

If  the  Taylor  series  of  a  function  f(x),  the 
general  term  of  which,  is 


(x-c)kfk(c) 

k*. 

is  truncated  after  n  terms,  the  resulting  expression  is 
a  polynomial  of  degree  <n.  Using  the  mean-value  theorem, 
which  states  that 


f  [ c  +  G (x-c ) ] 


f (x) -f (c) 

x-c 


for  some  0  satisfying  0<@<1,  we  now  consider 


(2.2.1) 


g(e)  =  pn(£)-(|E§-)n+1Fn(c) 


■ 


. 


■ 
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where 

(2.2.2)  F (?)  =  f(x)-fU)-  S 

n  k=l  K- 

Clearly  the  right  hand  side  of  (2.2.1.)  vanishes  at  |=x 
and  |=c .  Hence  by  Rolle's  theorem  its  derivative  with 
respect  to  £,  must  vanish  for  some  value  of  |  in  the  interval 
(x,c).  The  derivative  of  g(^)  is 


g'U)  = 


(n+l) (x-t)n 

(  \n+l 
(x-c) 


Therefore  we  have 


[ c+0 (x-c ) ] 


and  from  (2.2.2)  we  see  that  this  is  the  difference  between 
f(x)  and  the  Taylor  polynomial  and  is  equal  to 


V*>  =  [Sr|rlfn+1i  0+6 »• 

An  evaluation  of  the  Taylor  polynomials  yields 
the  following  points. 

1.  Within  the  range  of  convergence,  the  Taylor 

polynomial  yields  a  uniform  approximation  to  f(x) 
which,  improves  steadily  as  n  increases. 


13. 


2.  The  derivative  of  the  Taylor  polynomial  uniformly 
approaches  the  derivative  of  f(x)  as  n  increases. 

3.  Unless  we  can  eliminate  the  dependence  on  the 
derivative,  the  error  estimate  is  almost  unusable. 

4.  The  Taylor  polynomial  is  obviously  not  the  best 
approximation  to  f(x)  because  it  relies  on  the 
behavior  of  the  function  at  one  point  only. 

5.  The  necessary  data  for  the  Taylor  polynomial  of 
degree  n  consists  of  n+1  numerical  values :  f (x) 

and  its  first  n  derivatives  evaluated  at  the  single 
point  c . 

Prom  here  we  will  continue  on  to  much  more  accurate 
forms  of  approximation  -  those  that  use  information  from  more 
than  one  abscissa. 

Section  2.3  The  Lagrange  Polynomial 

The  ordinates  y^  involved  in  Lagrange's  form  of 
the  polynomial  L  (x)  =  y  (x)  of  degree  n,  which  takes 

on  the  same  values  as  a  given  function  f(x)  for  the  n+1 
distinct  abscissae  xq,x-^  .  .  .  ,xn,  are  displayed  explicitly. 

We  can  write  Ln(x)  directly  in  the  required  form. 


' 


. 
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(2.3.1)  Ln(x)  =  ^0(x)f(xQ)+. . .+in(x)f(xn) 

n 

=  2  i,(x)f(x,), 

k=0  K  K 

where  the  j}_k(x)  are  polynomials  of  degree  n,  to  be  determined 
by  the  requirement  that  L  (x^  =  f(x.)>  i  =  0,1,..., n.  Now 
the  expression  (2.3.l)  will  take  on  the  value  f(x^)  at 
x  =  x^  if 


i.  (x  .)  =  5 .  .  i  =  0,1 
i  J  lj 


i  n ,  j  —  0 , 1,.  .  . ,  n , 


where  6.  .  is  the  Kronecker  delta. 
lj 

Since  iL(x)  is  to  be  a  polynomial  of  degree  n 

which  vanishes  for  x  =  x  ,xn,...,x.  ,x.  ,x  ,  it  must 

o  1  l-l  1+1  n 

follow  that 


(2.3.2) 


n 

i .  (x)  =  c .  n  (x-x  .) 

1  j=o  ** 

jA 


where  c^  is  a  constant.  The  final  requirement  that  ^j_(x1)=l 
gives  us  the  following  expression  for  c. 


(2.3.3) 


c . 

1 


n 

n  (x. -x  . ) 

j=o  1  J 

jyi 


and  we  obtain  the  desired  influence  coefficients  ih(x)  by 


. 


J.  '  {  X  }  ,  X  $!~lV  fjW 


e  vi l  xbif;  do  i.riv’ 


) 
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introducing  (2.3.3)  into  (2.3.2). 

We  can  put  this  result  in  a  more  compact  form  by 
using  the  notation 


(2.3.4) 

n 

pn+1U)  =  n  (x-x  ). 
n  1  j=0  J 

The  derivative  of  pn+1(x)  is  expressible  as  the  sum  of  n+1 


terms  in 

each  of  which,  one  of  the  factors  of  P  ,-,(x)  is 

n+1 ' 

deleted . 

Thus  if  we  set  x  =  x^  in  this  expression,  we 

arrive  at 


Hence  we 

t  n 

p  , (x. )  =  n  (x.-x.)  =  l/c. 
n+1 v  l  y  '  l 

jyi 

can  write  (2.3.1)  in  the  form 

L  (x)  =  \  Pn+l(x)f(xk) 

k=°  (x-xk)P'+1(xk) 

^i(x)  is 

n 

It  is  obvious  that  the  coefficients  of  x  in 

l/Pl  (x. )  and  that  the  coefficient  of  xn  in 
n+1  n  L 

L  (x)  is 
n\  / 

th6n  +0f(x3>^'n+l(xi)- 

. 


. 


1 6. 


Section  2.4  The  Newton  Polynomials  and  Divided  Differences 


Let  L.(x)  denote  the  polynomial  of  degree  i 
which  passes  through  the  points  (xo,yo),  (x1,y1) ,  .  . 


Then  L  (x) 
n '  7 

can  be  written  as 

n 

(2.4.1)  Ln(x)  =  y0  +  .  (Li (x)  "  Ij1_1(x) )  • 

Now  by  definition  L. (x)-L.  n  (x)  vanishes  at  x  ,x-,,...,x.  n 
J  iv  l-l  o  1  l-l 

and  so  does 


(2.4.2) 

i-1 

n  (x-x  ). 

k=0 

Hence  the  L^(x) -L^_^(x)  and  the  product  can  differ  only  by 
a  constant  which  we  can  find  by  comparing  the  coefficients 
of  x1  in  the  product  and  in  L^(x)-L  -^(x).  The  coefficient 

of  x1  in  L . (x)  is 

l  ' 


(2.4.3) 

1  yk 

2  pi  rrr  r  -  y(x0>..*>xj)  • 

k=0  H+f  V  °  1 

Hence 


(2.4.4) 

Ln(x)  =  y  +  .ML  (x)-L  (x)] 

1=1 

yo  +(x-xo)y(xo,x1)  +  (x-xQ) (x-x1)y(xo,x1,x2) 

(2.4.4) 


■--!  '•  ■  :  1  .  .  wcV 


'  !£  '  ■'  '  '  a  ■  '  :  io 
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+  .  .  .  +  (x-xQ)  ( x— X-J^ )  .  .  .  (x-xn_]  )y(xo,x1,  . .  .  ,xn)  . 

Now  let  us  attempt  to  find  the  y(xQ,  .  .  „x^)  by  recursion. 
To  do  this  let  us  examine  the  linear  polynomial  defined  by 
(x  ,y  )  and  (x1,y1).  Using  (2. 4.4 )we  have 

(2.4.5)  L-l(x)  =  yo+(x-xo)y(xQ,x1) 

(2.4.6)  =  y1+(x-x1)y(x1,xQ) . 

If  we  put  x  =  x^  in  (£.4.5)  and  x  =  xq  in  (2.4.6)we  get 


yi'yO  ,  x  ,  x 

"x— 6T  = 

1  o 


Let  us  now  examine  the  quadratic  polynomial  defined  by 


(x1,y1) , (x2,y2)  and  (xQ,y0) . 


(2.4.7)  L0(x)  =  y1+(x-x1)y(x1,x2)+(x-x1)(x-x2)y(xo,x1,x2) . 


Putting  x  =  xq  in  (2.4.6)  and  (2.4.7)  gives  on  subtraction 


(xo-Xi)y(xo,Xi)  =  (xn-xn  )y(xn  ,xp)  +  (xn-xn  )(x0-x2)y(xo,x1,x2) 


o  lyj  v  1"  2'  '  o  l/x  O 


f 

V,-  i:  -- 


'  ■  ■  X  M  -) 
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so  that 


y(x0,x1,x2)  = 


y(x1,x2)-y(xQ,x1) 


x2  “  xo 


If  we  proceed  in  this  way  we  find  that 

y(x1,x2, . . . ,xi+1)-y(xo,x1, . . . ,x±) 


y(x0,x1, . . . ,x1+1 )  = 


i+1 


x .  .  ,  -  x 
1+1  o 


or  more  generally 
y(x1.x1+r....x1+k+1) 


y(x1+1,x1+2,  .  .  .;x1+lfJ.1)-y(xi.,xiJ.1„  .  .,x14.v) 


"i+k+1 


i*  i+1- 


"i+k- 


x .  . ,  .  -  x  . 

i+k+1  1 


These  expressions  are  called  "divided  differences" .  The 
polynomial  resulting  from  use  of  this  method  of  approximation 
is  the  same  polynomial  as  given  by  the  Lagrange  Method  and 
has  the  advantage  that  when  we  wish  to  find  an  approximating 
polynomial  of  higher  degree  we  need  only  extend  the  difference 
table,,  whereas  if  we  used  the  Lagrange  form  we  would  have  to 
modify  the  complete  system. 

Section  2.5  Error  Estimates  for  the  Lagrange-Newton 
Interpolation  Polynomials 

Let  us  emphasize  the  construction  of  the  approximation 
by  recalling  the  notation  yi=f(x^),  i=0,l, . . .,n.  The  polynomial 
constructed  from  the  ordinates  yields  at  best  an  approximation 
to  the  value  of  f(X)  at  any  point  X  not  equal  to  any  of  the  x^ . 


&&£  B>  9£>.T  ■;£  1:>  sr'-  >  rf  o  -  .  ,  w*.  .  ;Ic  f 

■  >£  m  1J:  S  .  ;■  ;  ;  -  ;t 

-  -  -  -  ■—-  ■  ■  -  . . . .  _  .  . 
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If  we  denote  the  value  of  f(X)  by  Y  and  the  value  of  L  (X) 
by  Yt  we  then  need  the  measure  of  the  deviation  Y-Y  . 

-Lj  l 

Let  Ln+1(x)  be  the  Lagrange  -  Newton  polynomial 
constructed  from  the  points  (xQ,y0) , (x  ,y  ),..., (xn>Yn) > 
(X,Y).  Then 

(2.5.1)  Ln+1(x)  =  yo+(x-xo)y(xQ,x1)+. . .+(x-xQ) (x-x1) . . . 

(x-xn_1)y(xQ. . .xn)+(x-xQ)(x-x1) . . . 
(x-xn)y(xQ,x1, . . . ,xn,X) , 

and 

Ln+l(x)  “  Y> 

and 

L  (X)  =  Yt . 
n  L 

If  we  subtract  (2 . 4 .4)  from  (2 . 5  •  l)  we  get 

Ln+i(X)  -  Ln(X)  =  (X-xo)(X-x1)...(X-xn)y(xo>x1,...,xn,X). 

However  this  is  not  of  much  use  in  its  present  form  because 
y(xQ,x^  , . . . ,xn,X)  uses  the  very  ordinate  Y  which  we  are 
trying  to  approximate.  Hence  let  us  try  to  put  the  error 


N  ii  ■■  1 V  r/ 


■ 

*~v  ••!  f.- 


■ 
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estimate  in  a  different  form.  To  do  this  we  make  the 
following  new  assumption;  let  the  (n+l)th  derivative  of 

f(x)  exist  in  some  interval  [a,b]  containing  the  and  X. 
Now  let  us  consider  f(x)  -  ^n+1(x)  which  vanishes  at  the 
points  x0,x1,...,xn,X.  Hence,  by  repeated  application  of 
Rolle’s  theorem,  the  first,  second,  ...  (n+1 )  th.  derivatives 
of  the  difference  vanish/  n+1,  n,...,l  times  respectively 
in  [a,b] . 


Let  |  be  the  point  where  the  (n+l)th  derivative 


vanishes . 


Then 


0. 


That  is 


(2.5.2)  fn+1(|)  =  (n+1 )  !  y(xQ,xi 


and  our  error  estimate  becomes 


Y-Y. 


(X-x0)(X-Xl)...(X-xn)fn+1(e) 


L 


(n+1 ) I 


This  leads  us  to  the  consideration  of  the  question 


what  is  the  best  polynomial  approximation?  -  which,  is 
considered  in  the  following  section. 


' 

as 

■  -  ■  r,  ,  J  ( x L  C  :a  >  .  .tsu  WoM 


X.  iqt  :  ;,X :  394 


' 
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Section  2.6  Chebyshev  Concepts 

Chebyshev  posed  the  following  question:  If  the 
degree,  n,  of  the  approximating  polynomial  is  fixed  in 
advance,  can  we  extract  from  the  set  of  all  polynomials 
of  degree  <n  that  polynomial,  p(x),  which  approximates 
f(x),  [f(x)  e  C[a,b]],  most  closely  in  the  sense  that 
max  |  f(x)-p(x)  |  is  least  in  [a,b]?  This  may  be  expressed 
in  a  more  mathematical  form.  If  p(x)  =  c  +c,x+...c  xn, 
then  the  greatest  of  the  deviations  |  f(x)-p(x)  |,  xe[a,b], 
is  a  positive  function  of  the  coefficients  c^,  which  we 
may  denote  by  p(c^).  If  we  allow  the  c^  to  vary  continuously 
over  all  real  values  then  p(c^)  has  a  least  value  and  we 
seek  the  polynomial  which  yields  the  least  value. 

We  shall  find  it  convenient  in  the  sequel  to 
abbreviate  this  statement  of  the  question  and  speak  of  the 
polynomial  of  approximation  p.  We  shall  denote  the  least 
value  of  p  by  minp .  Let  us  use  the  notation: 

ordinates  yQ,y1, . . . . ,yn+1 

abscissae  x  ,x,,....,x  with  the  ordering  x..  >x.. 

o  1  n+ 1  1+ 1  1 

To  fix  the  concept  let  us  take  a  special  case  of 
this  problem  first:  a  linear  polynomial  approximating 
three  points. 


■ 


x  l.bxsds 
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Let  the  ordinates  of  the  approximating  line  be 
yo“uo,yl’ul'  and  The  straight  line  with,  the 

ordinates  yQ-u0  and  y-^-u^  at  xq  and  x^  is 


(2.6.1) 


yM  = 


X-X. 

x  -x, (yo-uo}+ 
o  1  1 


x-x 

x, -lo(yl-ul} 


where  the  coefficients  of  the  ordinates  on  the  right  are 
the  Lagrangian  polynomials.  This  line  passes  through 

(x2'y2”U2^  and  hence 


y2-U2 


Xo-x,  x0-x 

(yo_uo)  +  (yi-ui} 


t>-  1 


o 


Now  multiplying  through,  by  (x^-xo)  we  have 


(2.6.2) 


(x2-x1) (y0-uQ) -(x2-x0) (y1-u1)+(x1-x0) (y2-u2) 


-  0 


or 


(2.6.3)  Uoqo"'ulcll+u2q2  =  J 


where  the  q^  are  the  positive  differences  of  the  abscissae  . 


It  is  convenient  to  express  the  u^  in  terms  of 

the  maximum  deviation  which,  occurs  at  one  or  more  of  the  x. 

1 

Thus  we  write  u.  =  pv.  where  |v. |<1  and  at  least  one  of  the 

1  r  q  I  I  *  — 


pfej.  c t:  d  a e ...  ■  ,  i  on  .  X  3 i  clT 


-  r;  '  X  ..X:  :b<  ■  f:  ..  .  "u-.i  -  3  id 
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v ^  =  +  1 .  We  can  regard  the  as  arbitrary  and  write 
(2.6.3)  as  an  equation  to  determine  p.  We  get 


P 


lyoqo~ylql+y2q2 

lVoqo-vlql+v2q2 


where  the  y^  and  q^  are  fixed  and  q^>0.  Only  the  v^  are 
variable  and  hence  are  at  our  disposal.  Clearly  the 
denominator  is  greatest  and  p  is  least  when  the  v^  alternate 
in  sign  and  each,  has  the  greatest  admissable  value,  unity. 
Hence 


min  p  = 


yoqo-ylql+y2q2 


Vql+q2 


This  shows  that  the  best  approximation  for  this  special 
case,  in  the  Chebyshev  sense,  is  the  one  which  has  the 
maximum  deviation  occurring  at  three  points  with  alternating 
sign. 


We  shall  now  consider  the  general  case  of  n+2 
ordinates,  y^, i=0, 1 , . . . ,n+l .  Let  us  denote  the  ordinates 
of  the  approximating  polynomial,  which  is  of  degree  n,  by 
y^-u^.  Now  the  polynomial  of  degree  n  which,  passes  through, 
the  n+1  points  (x^y^-u^)  ,i=0,l,  ...  ,n  is 


.  s  I  >  /.or 

a  ;  1  ■  '  ■ 

't:  w,,.a  Ms,: a 
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(2.6.4) 


Ln(x)  =  ^0(yi'UiUi(x) 


Where 

(2.6.5) 


n 


fn 


^  (x)  =  n  (x-x,  )  /  n  (x.-xj 
1  k=o  k  / k=o  1  k 


k^i 


k^i 


(the  Lagrangian  coefficients). 


Before  we  continue,  let  us  note  the  following 
two  points. 

1.  The  numerator  of  £.(x)  is  positive  for  x>x  . 

l  '  n 

n.  — ■  i 

2.  The  sign  of  the  denominator  is  (-1)  .  (In 

particular,  the  sign  of  ^j_(xn+i)  ls  (-l)n-1j. 

i=0, 1,  .  .  .  . ,  n,  that  is,  the  f.(x  .  )  alternate 
3  n+1 ' 

In  sign.) 

On  putting  x  =  xn+^  in  (2.6.4)  we  have 


(2.6.6) 


n 


^n+l~un+l 


=  2  (y .  -u  .  )i  .  (x  ,  ,  ) 

i=0  1  1  1  n+1 ' 


If  we  now  multiply  (2.6.6)  by  the  position  factor 


n 


q 


n+1 


n  (x.-x.) 

ij>  j=0  1  J 

i>j 


' 


■ 


' 
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we  can  write  (2.6.6)  in  the  form 
(2.6.7) 

n 

(yn+l-un+l)cin+l  =  il0(  -1) 

where  the  q^  are  the  product  of  all  the  absolute  value 
differences 

|xj~xk|  >  J>k>  j,k=0,l,  .  .  .  ,n+l. 

Writing  (2.6.7)  in  the  symmetric  form 


(2.6.8) 

n 

.  S^-l)1  n(y.-u.)q.  -  (y  ,  n -u  )  q  ,n  =  0 
i=0  l  i  i  wn+l  n+l/Mn+l 

and  u^  =  Pvj_>lv1|<  lj  we  have 


(2.6.9) 


y  0.  -y-i  cu  +  ...  +  ( -l)  y  q  +y  ,  q  n 
**0  0  Jl^l  v  '  Jn4n-Jn+1  m+1 

v  q  -V-.  qn  + .  .  .  +  (  -1  )nv  q  +v  ,  q  ~T~ 

ovd  1^1  v  '  nHn  n+1  vi+1 


It  is  now  obvious  that  the  minimum  value  of  p  occurs  when 

v.  =  1  and  the  signs  of  the  v.  alternate  since  the  y. 

1  i  1  °  l  J  i 

and  the  q^  are  fixed  and  are  not  at  our  disposal. 


We  can  draw  two  very  useful  conclusions  from 


I  i-r1 
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this  result.  Now  the  y.  and  q.  are  fixed  and  the  q. 

J 1  xi 

are  positive.  Hence  only  the  are  at  our  disposal  and 
|V^I<L  1*  Thus  we  have 


(2.6.10) 


ran  n  = 


yoqo-yl<1l+---  +  (-l)  yntt+A+l 
q0+qi+...+qn+qn+1 


The  second  conclusion  that  we  shall  draw  is  much,  more 
subtle.  Let  us  for  the  moment  go  back  to  the  u^  notation 
and  write  (2.6.8)  in  the  form 


(2.6.11) 


y  d  -yn  qn  +  .  .  .+y  . q  .  =  u  q  -unqn+.  .  .+u  .  n  q 

Jo^o  ^lk  -Jn+lHn+l  o  xo  1H1  -  n+lHn+l 


We  can  also  write 


(2.6.12) 


Now  the  preceding  analysis  shows  that  we  can  fit  a  polynomial 
of  degree  <n  through  the  n+2  points  (x^^y^-u^) ,i=0,l, . . . , 
n+1,  provided  only  the  u^  are  related  to  the  y^  by  (2.6.11) 
or,  what  is  equivalent,  provided  that  p  is  given  by  (2.6.9)- 
Hence  from  (2.6.12)  we  see  that  the  problem  of  approximating 
the  y^  by  a  polynomial  of  degree  <  n  is  equivalent  to  the 


: 


' 


.!  •>  I  ,:ir 


•  x  i  :  >  x } 


problem  of  approximating  the  u.  by  a  polynomial  of  degree 

n  . 


Let  us  denote  the  new  approximation  by  p1 
Having  done  so,  we  get 


(2.6.13) 


min  p! 


u  q  — u-,  q,  +...+( -1 )  u  q  ±u  .  ,q 
|  o^o  1^1  v  '  nHn  n+lHn+l 

q  +q-,  +  .  .  .  +q  +q  , 

^o  H1  Hn  Hn+1 


In  particular,  if  the  u^  alternate  in  sign  and  if  we  use 
the  notation  r.  =  |u.  we  have 

l  1  l 1 


(2.6.14) 


min  p1  = 


r  q  +rn  qn  + .  .  .  +r  nq  ,  , 
o^-q  InL _ n+1  Hn+1 

q  tqn  + .  .  .  +q 
Ho  ^1  Hn+1 


Since  the  problems  are  equivalent  we  have 

min  p1  =  min  p. 

From  (2.6. l4)  we  see  that,  in  effect,  min  p 5  lies  between 
the  greatest  and  the  least  values  of  r^,  that  is,  it  is 
a  weighted  average. 

These  two  conclusions  lead  us  to  the  following 


theorems . 


(c 
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Theorem  (2.6.1)  There  exists  a  unique  polynomial  of 
degree  £n  which  furnishes  the  closest  approximation 
to  a  set  of  n+2  distinct  ordinates.  The  deviations 
between  the  ordinates  of  the  approximating  polynomial 
and  the  prescribed  ordinates  alternate  in  sign  and  are 
of  equal  magnitude . 

Theorem  (2.6.2)  If  a  polynomial  of  degree  <  n  furnishes 
deviations  u^  (from  n+2  prescribed  ordinates)  of  alter¬ 
nating  signs,,  but  unequal  magnitude,  then  the  closest 
approximation  to  the  prescribed  ordinates  lies  between 
the  greatest  and  the  least  of  the  absolute  value  of  the 
u^  and  is  the  weighted  average  of  the  u^;  the  weights 
being  the  q^ . 

We  can  now  return  to  the  original  problem  by 
generalizing  theorem  (2.6.1)  in  two  steps.  First  we 
consider  the  best  approximation  to  a  set  of  m  ordinates 
(m>n+2) .  Then  we  consider  the  case  in  which  the  ordinates 
may  have  any  abscissae  in  a  continuous  strip  [a,b].  The 
preceding  work  enables  us  to  formulate  the  following  theorem. 

Theorem  (2.6.3)  The  best  polynomial  approximation  of 

degree  <  n  to  a  finite  set  of  m(m>n+2)  points  is  that 

in 

one  corresponding  to  the  greatest  of  the  (n+p)  values  of 
min  p  at  n+2  points  selected  arbitrarily  from  m  points. 
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Proof : 

Let  us  use  the  notation 

a  for  the  greatest  of  the  values  of  min  p, 

E  for  the  set  of  n+2  points  which  furnishes  o, 
y^  for  the  set  of  ordinates  of  these  points. 

If  the  absolute  value  of  the  deviation  u,  at  any  point  p 
(of  the  set  of  points)  not  belonging  to  E,  is  greater  than 
o,  then  let  E1  be  the  set  of  points  obtained  by  deleting 
one  of  the  points  of  E  and  adding  p,  the  deletion  being 
made  in  such  a  way  that  deviations  of  the  points  of  E 1 
alternate  in  sign.  Write  Up  =  a+ i,  t>0 .  Then  from 
theorem  2 

min  o  =  (c+-T)qyaqi+...+a!jI>+1 

Ft  ®  > 

*  qT+q’+.  .  .+q!  ,  , 

o  1  n+1 

which  is  a  contradiction  that  proves  the  theorem. 

We  also  have  the  following  corollary  to  theorem 

(2.6.3) . 

Corollary :  If  f(x)  is  continuous  in  [a*b],  the  closest 

polynomial ' approximation  of  degree  <  n  is  that  which 
corresponds  to  the  greatest  of  the  values  of  min  p  on 
any  set  of  n+2  points  of  [a^b] . 

Proof : 


The  best  approximation  over  a  set  of  n+2  points 


t 

, 
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in  [a,b]  is  given  by  (2.6.10).  This  formula  demonstrates 
that  p  is  a  continuous  function  of  the  x ^  in  the  set  in 
[a,b] .  Now  this  function  admits  a  maximum  p  and  attains 
it  for  a  particular  set  E  in  [a,b] .  Hence  we  can  prove 
the  corollary  by  the  same  argument  as  that  employed  to 
prove  theorem  (2.6.3)- 
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CHAPTER  III 

RATIONAL  APPROXIMATION 


Section  3.1  Introduction 


In  the  last  chapter  we  observed  that  there  is 
always  a  polynomial  approximation ,  of  degree  <  n,  of  best 
fit  (in  the  Chebyshev  sense)  to  a  set  of  n  +  2  points. 

In  this  chapter  we  will  investigate  rational  approximations 
to  continuous  functions.  We  will  then  proceed  to  the 
difficult  problem:  Using  a  similar  line  of  attack  to  that 
of  Chapter  II  section  6  can  we  find  a  min  p  for  a  rational 
approximation  of  the  form  p  (x)/q~(x),  where  pn(x)  is  a 
polynomial  of  degree  <  n  and  qm(x)  is  a  polynomial  of 
degree  <  m^  to  a  given  set  of  n  +  m  +  2  points?  Also, 
does  a  best  rational  approximation  always  exist? 


Section  3.2  Rational  Approximations  to  Continuous  Functions: 
Statement  of  the  Problem 

Chebyshev ls  formulation  of  a  class  of  problems 
of  this  nature  is  admirable  and  lucid  and  we  follow  him 
here  . 


It  is  required  to  determine  a  function 
involving  n  +  1  arbitrary  parameters  p^,  P2>  •  • 

such  that  the  greatest  deviation  |  y(x)  -  f(x) 
y (x)  and  the  given  function  f(x)  in  the  interval 
-  a  <  x  <  a  is  least.  More  precisely,  we  require 

min  max  I  y(x)  -  f(x)  I 
pi 


y(x) 

,3  Pn+1 
between 


-a  <  x  <  a 


“-a  1o 


. 
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This  formulation  embraces  the  following  problems 


1.  To  minimize  with  respect  to  the  Pj_ 

(3.2.1) 

r(x)  =  |  pQxn'v  +  p1xn'v_1  +  ...  +  pn_v  -  f(x)  | 


in  -  a  <  x  <  a 


2.  To  minimize  with  respect  to  the  p^ 

(3.2.2) 

/  \  ,Po 

r(x)  = 


n-v  ,  n-v-1  . 

p0x  +  P]X  +  . . 


•  +  Pn- 


m-Li  .  m-^L-1  . 

q  x  ^  +  q.x  ^  +  .  . 

Lo  ^1 


in  -  a  <  x  <  a.  The  qn-  are  given. 


n-  v 


•  +  q 


m-q 


-  f(x) 


3.  To  minimize  with  respect  to  the  p.*  and  q. 


(3.2.3) 


r  (x)  = 


P0X 


n-v  .  n-v-1 

+  PqX 


•  •  +  Pn- 


n-  v 


q0xm“P  +  q^x™"^-!  + 


+  ^m-ii 


-  f(x) 


in  -  a  <  x  <  a . 


This  is  sufficient  generality  for  our  present 
purpose.  Problem  1  is  a  special  case  of  problems  2  and 
3.  Problem  2  has  a  fixed  denominator  in  the  rational 
fraction  and  we  shall  not  consider  it  further.  Problem  3 
embraces  problem  1  and  problem  2.  Here  we  shall  consider 
only  problem  3. 


One  of  the  fundamental  theorems  concerning 


qa~n  Aw  l<  ■  :  mi'.n.i  r 


33. 


problem  3  is  usually  quoted  In  the  form  given  by  Achiezer; 
however,  a  different  statement  of  the  theorem  was  given 
by  Chebyshev,  which  is  In  some  respects  more  lucid.  More¬ 
over,  Chebyshev1 s  proof  of  the  theorem  provides  a  good 
understanding  of  the  theory. 

We  will  consider  first  Achiezer 's  form  of  the 
theorem  and  then  Chebyshev' s.  The  two  forms  have  differ¬ 
ent  contents  -  in  Achiezer *s,  the  number  of  extrema  for 
a  polynomial  pair  of  orders  n  -  v,  m  -  q  is  defined  to  be 
>  m  +  n  +  2  -  min  (v,p),  while  in  Chebyshev* s  the  number 
of  extrema  is  fixed  and  the  theorem  gives  information 
about  the  greatest  possible  orders  of  the  polynomial  pair. 

Section  3-3  Achiezer* s  Form  of  Chebyshev  *  s  Theorem 

We  will  first  give  a  number  of  theorems  prepara¬ 
tory  to  theorem  3-3-2  which  is  our  final  result. 

Section  3-3.1  A  Generalization  of  de  la  Vallee  Poussin's 
Theorem 

Before  stating  the  generalization  let  us  recall 
to  mind  de  la  Vallee  Poussin's  theorem  which  may  be  stated 
thus  . 


Let  y±,  i  =  0,  1,  . . . , 


n+1,  be  a  set  of  n+2 


>  ;  ■  .  i ;  ;  . :  .  cm' 
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ordinates  at  the  abscissae  xq  and  let  y(x)  be  a  polynomial 
of  degree  <  n.  Let 

r±  =  y(xi)  -  y±  1  =  0,  1,  .  .  . ,  n+1 


and  let  the  r ^  alternate  in  sign.  Then  the  polynomial  of 
best  approximation  (of  order  <  n)  on  this  set  of  points 
has  a  deviation  which  lies  between  the  greatest  and  least 
of  the  absolute  values  of  the  r . . 


The  significant  point  here  is  that  the  r^  must 
alternate  in  sign.  The  lower  bound  is  useful,  the  upper 
bound  less  so  since  in  approximating  to  a  continuous 
function  we  must  consider  all  possible  sets  of  n  +  2  points. 

With  this  preamble  we  proceed  to  the  generalization. 
We  use  the  notation 


(3. 3. l.l) 

Q(x) 


n- v  n- V- 1 

PQx  +  Pjx  +  . . .  +  pn_ y 

q0xm_^  +  q-]_xm“1J'“-1'  +  .  .  .  +  qm_|JL 


=  s(x) 


p(x) 

q(x) 


Here  s(x)  is  a  continuous  function  which  never 
vanishes  in  the  closed  interval  [-a, a]  ,  m  and  n  are 
given,  and  0<ir£.m,  0<v<n  .  The  point  of  including 

the  parameters  \i  and  v  is  that  we  wish  to  consider  all 
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polynomials  p(x)  of  degree  <  n  and  all  polynomials  q(x) 
of  degree  <  m.  No  special  notation  is  needed  for  this 
since  the  coefficients  of  p(x)  and  q(x)  -  in  particular, 
the  leading  coefficients  -  may  have  any  finite  values, 
zero  or  non-zero.  It  is,  however,  convenient  in  the 
analysis  which  follows  to  indicate  explicitly  the  degrees 
of  a  particular  p(x)  and  particular  q(x)  by  using  the 
parameters  \i  and  v. 

We  assume  that  p(x)  and  q(x)  have  no  common 
divisor  and  that  q(x)  does  not  vanish  in  [-a, a] . 

Let  f(x)  be  continuous  in  [-a, a]  and  let  the 
maximum  absolute  deviation  between  f(x)  and  Q(x)  in  [-a, a] 
be  that  is, 

pQ  =  max  |  Q(x)  -  f(x)  | 

-a  <  x  <  a 

We  can  now  state 

Theorem  3. 3. 1.1:  A  generalization  of  de  la  Vallee 
Poussin 3 s  theorem. 

Let  Q(x)  =  s(x)  p- S-x- 1  remain  finite  in 

qix ) 

[-a, a]  and  let  the  deviation 

Q(x±)  -  f(x±) 

at  the  consecutive  points  x-^,  x^,,  ...  ^  x^  of  [-a, a] 
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assume  the  values 

v  -  V  '  f-1)”'1  \ 

where  the  A^  are  positive  and  where  N  =  m  +  n  +  2-  min(p,v)  . 
Then  p  >  min  A^  i  =  1,  2^  ...  N. 


Proof 


Let  R(x) 


be  such  that 


s(x) 


P0*  xn  +  Pl>  xn_1  +  ...  +  pn! 
qo>  xm  +  qi»  x111’1  +  ...  +  qm' 


pR  <  min  A^  . 


Then  we  can  prove  a  contradiction.  For  writing. 


A(x)  =  R(x)  -  Q(x)  =  - [  Q(x)  -  f(x)]  +  [R(x)  -  f(x)] 


we  see  that  the  sequence  of  numbers 


A(xx) ,  A(x2) ,  . . .  ,  A(xn) 

are  different  from  zero  and  alternate  in  sign.  Since  A(x) 
is  continuous  in  [-a, a],  it  follows  that  A(x)  has  at  least 

N-l  =  m+  n+  l-  d  (d  =  min(p,v)) 

zeros  in  [-a, a].  But  that  is  impossible  since 

A(x)  =  s(x)  p8(x)  q(x)  -  p(x)  q'(x) 

q(x)  q 1 (x) 


1 
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and  the  degree  of  the  numerator  is  at  most  m  +  n  -  d 
(at  most,  since  we  admit  that  the  leading  coefficients 
of  p'(x)  and  q'(x)  may  be  zero).  This  contradiction 
proves  our  theorem. 

Corollary  to  Theorem  3.3.1.!.  At  this  point  we  recognize 
that  we  cannot  do  better  if  we  can  fi'nd  a  Q(x)  whose 
deviation  from  f (x)  at  a  set  of  points  x^,  i  =  1,  2,  . . .  N, 

are  A,  -A,  ...  ,  (-l)N_1  A.  For  Theorem  3-3. 1.1  assures 
us  that  no  rational  function  of  the  form  Q(x)  can  have 
a  smaller  maximum  deviation.  For  practical  purposes  this 
is  sufficient;  but  we  can  go  further  and  show  that  the 
Q(x)  of  minimum  deviation  is  unique. 

We  are  assuming  existence  here.  For  completeness 

we  state: 

Theorem  3- 3. 1.2:  The  Existence  Theorem.  Among  the  functions 
Q(x)  there  exists  at  least  one  function  for  which  is  a 
minimum . 


However  Q(x)  includes  the  polynomials  of  degree 
n-v.  Hence  the  existe.nce  of  best  rational  approximations 
with  a  polynomial  denomination  is  not  assured  by  the  above. 

We  now  proceed  to  the  final  theorem  of  this 
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Section  3-3.2  Achiezer’s  Form  of  Chebyshev 1 s  Uniqueness 

Theorem 


The  function  Q(x)  of  the  form 


Q(x)  =  s ( x ) 


n-v  .  n-v-1  . 

PQX  +  p-jX  +  ...  +  p 


n-  v 


q  xm^  +  qn xm“^-1  +  ...  +  q 
Ho  Hm-q 


which  deviates  less  in  [-a, a]  from  a  given  continuous 
function  f(x)  than  any  other  function  of  the  same  form 
is  completely  characterized  by  the  following  property. 

The  difference  Q(x)  -  f(x)  assumes  the  maximum  values 
±  Pq  not  less  than  m  +  n  +  2  -  d  times  in  the  interval 
[-a, a],  positive  and  negative  deviations  occurring  alter¬ 
natively. 


In  the  statement  of  this  theorem,  we  have 
supposed  that  p(x)/q(x)  is  irreducible  -  that  is,  no 
linear  factors  of  the  numerator  correspond  to  the  linear 
factors  of  the  denominator.  The  proof  of  this  theorem 
is  quite  delicate.  It  seems  simpler  to  prove  the  theorem 
first  for  a  polynomial  and  then  modify  the  proof  to  cover 
the  rational  function  Q(x) .  The  proof  for  the  polynomial 
case  is  due  to  de  la  Vallee  Poussin. 

The  theorem  for  a  polynomial  p(x)  of  degree  <  n 
amounts  to  the  assertion  that  the  minimum  polynomial  p(x) 
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has  the  property  that  r(x)  =  p(x)  -  f(x)  achieves  its 
greatest  value,  -  p,  at  least  n  +  2  times  in  the  interval 
[-a, a],  positive  and  negative  deviations  occurring 
alternately.  Let  x^,  x^,  ...  x^  be  the  abscissae  of  the 

points  at  which  r(x)  =  -  p .  We  prove  the  theorem  by 

showing  that  N  <  n  +  1  involves  a  contradiction.  Since 
r(x)  is  continuous,  we  can  find  a  small  interval  6^ 
enclosing  each  x^  such  that  r(x)  is  one-signed  in  each 
interval  5^. 

First  let  us  notice  that  if  all  of  the  maximum 
deviations  are  one-signed,  we  can  achieve  a  better 
approximation  by  adding  to  p(x)  a  small  constant  of  the 
same  sign  as  -r(x) . 

Suppose  then  that  the  sequence  of  maximum 
deviations  comports  k  <  n  groups  of  consecutive  terms  of 
contrary  sign.  Let  a  change  of  sign  occur  between  the 
intervals  6.  and  6..,.  Clearly  6.  and  <5 .  .  are  not 

1  1+ 1  1  It  1 

contiguous  (since  r(x)  does  not  vanish  in  either)  and  we 
may  choose  a  point  between  6^  and  with  the  terminal 

points  of  the  intervals  being  excluded.  (See  figure  1) 


8  -  -  :  '  ■  ,  ';,v  :r. 
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Graph  of  r(x)  =  p(x)  -  f(x) 

Figure  1 

Now  let  ...  be  the  set  of  points 

so  chosen  corresponding  to  the  k  consecutive  groups.  Let 

k 

&W  =  [-Sgn(r(x  ))]  n  (£.-x)  , 

1  i=I  1 

This  follows  from  the  fact  that  we  want  sgn[p'(xi)  ]  = 
-sgn[r(x^)].  Then  the  polynomial  ^(x)  has  the  same  sign 
as  -r(x)  in  each  of  the  intervals  6^  and  does  not  vanish 
there.  Since  k  <  n.,  0(x)  is  a  polynomial  of  degree  <  n  and 

p(x)  +  e^(x)  (e  small  and  positive) 

is  a  polynomial  of  degree  <  n  which  deviates  less  from 
f (x)  in  each  interval  6^  than  does  p(x) .  Now.,  if  we 
can  show  that  it  deviates  less  in  the  remainder  of  the 
interval  [-a, a],  we  shall  have  proved  our  contradiction. 

To  complete  the  proof:  let  p *  <  p  be  the 
bound  of  |  r(x)  |  in  the  intervals  of  [-a,a]j  exterior 
to  the  intervals  6^. 


Let  e  be  chosen  so  that 


■  \:C'. ...  0‘;  xo;0,i:  . ,  :b^o 


4l . 


I  |  <  p  -  p '  in  [-a, a]  ; 

then  |  f(x)  -  p(x)  -  e^(x)  |  <  |  r(x)  |  +  |  ep'(x)  |  and 
for  x  not  included  in  6^,  i  =  I,  II.,  .  ..,  k  we  have 
I  f(x)  -  p(x)  -  ep'(x)  |  <  p{  +  (p-pl)  =  p 

or  |  f (x)  -  p(x)  -  €0(x)  |  <  p  . 

This  is  a  contradiction  because  we  assumed  that  p(x)  was 
the  best  approximation  with  p  =  min  p. 

The  modification  of  the  proof  to  deal  with  a 
rational  function  Q(x)  requires  some  minor  changes.  We 
now  achieve  a  contradiction  by  assuming  that  the  number 
of  maximum  deviations  of  r(x)  =  Q(x)  -  f(x)  is 
N<m+n+l-d  (d=  min  (q,v)).  The  polynomial  $  {x) 
is  chosen  as  before  but  its  degree,  k,  is  now  determined 
by  k  <  N  -  1;  that  is,  k  <  m  +  n  -  d.  We  now  choose 
polynomials  a(x)  and  b(x)  such  that 

0{x)  =  q(x)  a(x)  +  p(x)  b(x), 

where  the  polynomial  a(x)  is  of  degree  <  n  and  b(x)  is 
of  degree  <  m. 

Now  let  us  consider  the  comparison  function 
Q’(x),  instead  of  Q(x) .  We  have 


Q’(x)  =  Q(x)  +  s(x)  — j — r-, — -  - , — pp 
v  v  v  ’  q(x) [ q(x) -  eb(x)J 


•  'i: 

asdraii/n  -fi'  |snl 

.  (  V  t  ,.\  I tt'il  ■■  b  )  ' 

'  : '  ■  ■  2:  i:  -  r,  :  -  If.  ■ '  ,  :  . .ir  el 
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The  reasoning  applied  in  the  polynomial  case 
goes  through  unchanged  if  we  add  the  extra  condition  that 
e  be  chosen  such  that: 

q(x)  -  eb(x) 

has  the  same  sign  as  q(x)  in  [-a, a] . 

The  only  point  remaining  to  be  established  is 
that  04 (x)  is  of  the  same  form  as  Q(x) .  But 


Q'(x) 


e[q(x)a(x)  +  p(x)b(x) j  4 

q(x) [ q(x)  -  €b(x)j  / 


=  s(x) 


{ 


P ( x ) [ q ( x ) — eb ( x ) ]  +  e[q(x)a(x)+p(x)b(x)  4 

q(x) [ q(x) - eb (x) ] 


S(x)  1 

lq(x)  -  eb(x)  J 


which  is  of  the  same  form  as  Q(x) .  Thus  we  reach  the 
conclusion  that  Q*(x)  has  a  smaller  deviation  than  Q(x) 
which  contradicts  our  original  assumption. 


Section  3.4  Chebyshev * s  Theorem 

We  will  use  the  notation  of  (3.3.1.1)*  with 
qm  =  1,  for  Q(x)  in  what  follows. 


:  :  O  '  !  '  d  5  (  d  '  »n :  -T  ■  :.C  .■  ;  u!  IT 


.  -D..; 


■  '  '  i  3  a  L  ; :  OX'  !w 


..  v  :  ■  .  <  ...  -v;-  J  :J. ,,  dd 


r  .'•  .  •  .  'v  ' "  ID  ' ' .  '  •  ?o.\;  :)€■: 

.  "  -  — 
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Chebyshev's  theorem  can  be  stated  as  follows: 

If  the  difference 

Q. (x)  -  f (x)  =  r(x) 

achieves  its  greatest  absolute  value  L  less  than  m  +  n  +  2 
times  in  [-a, a].,  then 

p0  =  Pi  =  ...  =  pd_i  =  0 

qG  =  Qi  =  •••  =  qd-l  =  0 

where  the  number  of  extrema  in  [-a, a]  is  m  +  n  +  2  -  d. 

It  is  important  to  note  that  part  of  Chebyshev's 
analysis  is  restricted  by  the  assumption  that  f(x)  (the 
function  to  be  approximated)  is  differentiable.  His 
reasoning  was  -  in  part  -  based  on  the  fact  that  for  a 
differentiable  f(x),  the  extrema  are  the  common  zeros  of 

r2(x)  -  L2  =  0 

and  (x+a) (x-a)  r'(x)  =  0 

in  [ -a, a]  . 


For  f(x)  dif ferentiable,  there  is  no  theoretical 
difficulty  when  the  number  of  extrema  is  m  +  n  +  2  .  For 
the  two  equations  above  imply  at  least  m  +  n  +  2  equations 
for  the  determination  of  the  n  +  1  quantities  p^,  i  =  0, 

.  .  . n,  and  the  m  quantities  q^*  i  =  0,  ...  ,  m-1,  and 
the  deviation  L.  There  are  in  principle  enough  equations 
for  the  solution  of  the  problem. 


.Jb  -  ■£  +rn  •+  ff!  8 .!:  ;BtB- 


. 


L  .  i  ..  •  "  '  ■  .  ■  '  '  .!'/■  i:lb 
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Let  us  now  consider  the  case  where  the  number 
of  extrema  is  less  than  m  +  n  +  2.  Although  this  case 
is  more  difficult,  it  is  possible  to  show  that  it  is 
soluble  in  principle.  Supposing  that  the  number  of 
extrema  is  a  (<  m  +  n  +  2),  we  shall  presently  set  up  a 
system  of  m  +  n  +  1  equations.  When  a  solution  exists, 
it  is  possible  to  obtain  from  these  -  by  a  process  of 
elimination  -  another  set  of  m  +  n  +  2  -  a  equations  for 
the  determination  of  the  m  +  n  +  2  +  a  quantities  p^,  q^, 
L,  x  .,  j  =  1,  . . .  ,  a.  Putting  x  =  x  ^  in  both  of  the 
equations 

r2(x)  -  L2  =  0 

(x+a) (x-a)  r 1  (x)  =  0 

we  have  2a  more  equations,  ie.  m+n+2+ain  all, 
which  is  sufficient. 


Proof  of  the  Theorem 


In  the  proof,  we  shall  consider  only  the  case 
a  <  m  +  n  +  1.  We  will  use  the  following  notation: 


Q(x)  = 


pQxn  +  P]_xn~ 1  +  .  .  .  +  pn _ 

qQxm  +  q]_xm  1  +  ...  +  qm_]_  x  +  1 

n 

.2 


=  i=o  Pn-i  x' 


m-1 

1  +  .2  q„  n  _.  x 


i+1 


1=0 


^-.m-l-i 


J 


' :  '  '■  . 
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x^,  i  =  1,  ...  ,  a,  for  the  extrema  of 


and 


r(x)  =  Q(x)  -  f ( x ) 

L  =  .max  |  r(x)  | 

-  a  <  x  <  a 


We  need  a  number  of  preliminary  lemmas. 


Lemma  1 .  The  quantity  L  is  not  reduced  to  its  smallest 
value  if  the  system  of  equations 


(3.4.1) 


a 

Z 

j=l 


ij 


A 


0 


where 


9r(Xj) 

SPi-1 

dr  ( x  n- ) 
^i-n-2 


i  =  1,  ,  n+1 

i  =  n+2,  ...  ,  n+m+1. 


has  no  solution  other  than  A  . 

J 


0,  j  =  1, 


a . 


The  proof  of  Lemma  1  requires  two  further 
lemmas.,  1.1  and  1.2  below.  However,  let  us  note  in 
passing  that  if  a  solution  exists  in  which  one  or  more 
of  the  Aj_  are  non- zero,  then  we  can  use  a  of  the 
equations  to  yield  a  determinantal  equation  for  the 
p^,  q^,  x^,  and  L.  This  equation,  together  with  the 
remaining  m  +  n  +  1  -  a  equations  constitutes  a  system  of 


. 

- 


* 


'•  J  (xV'i  j 


- 
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(n  +  m+  l-a)  +  l  =  n  +  m+  2-  a 


equations  as  noted  above. 


Lemma  1.1.  If  the  system  (3 .4 .  ])  has  no  solution  other 

than  X.  =  0,  j  =  1,  ...  ,  o ,  then  we  can  find  a  set  of  numbers 

J 

r|^,  i  =  !_,  .  .  . ,  n+m+1 *  such  that  the  a  equations 


(3.4.2) 


n+m+1 

ill 


f 


ji  + 


r(xj) 


are  satisfied 


Proof  of  Lemma  1.1 

There  are  two  cases  to  consider: 

1)  a  =  m  +  n  +  1.  Then  the  matrix  of  the  f. .  is 

ij 

square  and  of  rank  m  +  n  +  1.  Hence  the  equations 
(3 . 4 . 2)  certainly  have  a  non-zero  solution. 

2)  o  <  m  +  n  +  1.  Then  the  matrix  of  the  f .  .  is  of 

ij 

rank  a.  We  may  then  attribute  arbitrary  values  to 
(n+m+1)  -  a  of  the  r\^  and  solve  (3.4.2)  for  the 
remaining 


Lemma  1.2.  There  exists  an  e  >  0  such  that 

max  |  rQ(x)  |  <  L  , 

-a  <  x  <  a 


rQ(x)  =  Qq(x)  -  f(x) 


where 


•1 


.Ls  o:i  r. 


■ 


. .  • 
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and  Q  (x)  is  the  rational  function  which  is  obtained 
o 

when  the  coefficients  p^,  q^,  in  Q(x)  are  replaced  by 

Pi  -  eTi±,  q±  -  en±  . 


Proof.  We  can  divide  [-a,a]  into  two  sets  of  intervals, 

s^  and  s2,  where  s^  consists  of  intervals  enclosing  the 

extrema  such  that  r(x)  does  not  vanish  in  any  of  these 

intervals,  and  s2  consists  of  the  remaining  intervals. 

Let  us  consider  a  typical  interval  of  s^.  At 

the  extremum  x  .,  we  have 
J 


r 

o 


rUj) 


m+n+1 
€  2 
1=1 


ar(xj) 

‘l  -N 

3P±-1 


+  0(e2) , 


(where  we  are  assuming  p  .  =  q.),  since  Q(x)  -  and 

n_t_  J  J 

hence  r(x)  -  is  a  differentiable  function  of  the  parameters 
p^  so  long  as  the  denominator  of  Q(x)  does  not  vanish  in 
[-a, a].  It  is  easy,  but  messy,  to  prove  that  the  term 
0(e)  can  be  handled.  Hence  we  shall  ignore  it.  By 
reason  of  Lemma  1.1  we  have 


rQ(xj)  =  r (x j)  [  1-e]  +  0(e2)  <  r(xj.) 
since  e  >  0. 


Now  let  x  be  any  other  point  of  this  interval 


a  lo  lav  o:;  Jao;-'  r  n  -i 
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and  let  the  Interval  be  chosen  so  that 

r(x)  =  r(xj)  +  61 

which  can  be  written  as 

2  n  Mill  =  2  T).  +  5 

1  3p±  1  1  3  Pi  2 

since  r(x)  and  dr(x)/dp^  are  continuous  functions  of  x. 

Now 

rQ(x)  =  r(x)  -e[2T]1  +  0(e2) 

=  r(x)  -€[2t)±  +  62]  +  0(e2), 

that  is, 

rQ(x)  =  r(x)  -  e[r(xj.)  +  62]  +  0(e2) 

=  r(x)  -  e[r(x)+61+62]  +  0(e2) 

=  r(x)[l-e]  -  £(6^62)  +  0(e2) 

We  require  that  |  6^  +  62  |  <  1  so  that  we  have 
|  rQ(x)  |  <  |  r(x)  |  and  therefore  |  rQ(x)|  <  L. 

We  have  yet  to  fix  e.  If  we  let  |  |  r(x)  |  -  L  |  <  5 
in  any  interval  of  s2,  then  it  is  sufficient  to  choose  e 
such  that  |  rQ(x)  -  r(x)  |  <  6  for  all  x  contained  in  s2 . 


■ 
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This  completes  the  proof  of  Lemma  1.2  and  hence  the  proof 
of  Lemma  1 . 


We  are  now  in  a  position  to  prove  Chebyshev  *  s 
theorem.  We  will  first  set  up  the  specific  form  of 
equations  (3 .4 . 1)  when  r(x)  =  Q(x)  -  f(x)  and 

,  .  =  Po*n  +  Pl^'1  +  •  •  •  +  Pn 

q0xm  +  q^x”-1  +  .  . .  +  qm_1x  +  1 
n  n-i 

Jo  Ppx 


m-1  m-i 

1  +  ilo  qiX 


where  the  abscissae  of  the  extrema  of  r(x)  are  x  • 
j  =  1,  2,  .  .  .,  a,  where  a  <  m  +  n  +  2.  We  have 


dr(x) 

dp. 


n-i 


n-i 


x 


q(x)  q  (x) 


1  =  0,1,  .  .  . ,  n. 


and 


dr  (x) 
dq. 


m-i  /  \ 
x  P(x) 

q2(x) 


i  =  0,  1,  ...  ,  m-1 


Hence  the  equations  (3 . 4  .  l)  become,  in  this  case. 

At 


a 

Z 


xk 


k=l  q2 
a  A 


qHxk) 


q(xk)  xk 


n-  i 


=  0,  i  =  0,  1,  .  .  . ,  n 


k 


k=l  J2 


q  (xk) 


p(xk)  xk 


m-i 


=  0,  i  =  0,  1, 


m-1 


i  arfe^GfarlO  evo^rq  op  p®l#lLoq  ■■  '  won  ■  -t.b  o  . 
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Now,,  if  we  multiply  the  first  set  of  equations  by  p  , 
i  =  0,  1,  . .  n,  and  the  second  set  of  equations  by 
i  =  0,  1,  .  ..,  m-1,  and  subtract  we  have 


(3.4.3) 


l  _3- 


“  -Z? - \  ¥(xj  =  0 

k=i  q  (%)  k 


n 


m-1 


where  ¥(x)  =  q(x)  ±ZQ  Pn_±x  -  p(x)  ±ZQ  Ym_1_1: 


.i+1 


It  is  most  convenient  to  set  out  the  proof  as 
a  number  of  lemmas . 


Lemma  2 .  If  Y(x)  is  of  order  a  -  1  and  if  the  coefficients 
^iJ  ^i  can  c^osen  so  that 

(3.4.4)  T(x)  =  (x-x2) (x-x3) . . . (x-xg) 

a 

-  n  (x-x . )  , 

j=2  J 

then  A^  =  0 . 

Proof .  The  lemma  is  obvious  from  the  form  of  (3-4.3). 

If  T(x)  is  of  the  form  (3-4.4),  then  ¥(x^)  vanishes  for 
k  /  1.  It  follows  from  (3-4.4)  that  A-^  =  0. 

A  similar  lemma  would  enable  us  to  show  that 

A^  =  0  if  Y(x)  is  of  the  form  ^(x),  that  is,  a  product 

of  all  the  factors  (x-x.),  j  =  1,  2,  ...,  a,  j  ^  k. 

J 


.  G  - 
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Lemma  3 .  ¥(x)  can  be  reduced  to  any  of  the  forms 

¥^(x),  k  =  1,  2,  a  as  long  as  one  or  more  of  the 

coefficients 


Pp*  • • •  j  Pd_x 


d  =  n  +  m  +  2-  a 


q 


o ’  M1J  * ’ 


J  •  •  •  J 


q 


d- 1 


is  non- zero. 


Proof  .  We  may  write 


.  .  .  n  -i  m—  1 

(3-^.5)  ¥(x)  =  q(x)  .2  (3  .xx  -  p(x)  .2  y  ,  .x 

w  v  '  41  ;  i=o  n-i  '  i=o  rm-l-i 


i+1 


However, 


q(x)  .2  (3  . (x)  = 

'  i=o  n-iv  ' 


m-1 

j=o  ^m-j-r 


n 


[1  +  .2  q  .  nxJ*+1]  ,2„  ^  ,-x1 


i=o  n-i 


n 


m+n-1  min (^, m-1) 


=  .2  p  .x  +  „ 2  .2 

i=o  n-i  l=o  j=o 


^n-l+  j^m-  j-lJ 


1+1 


and 


m-1 

p  ( x )  .2  Y  ..  .xJ+1 

'  j=o  'm-l-j 


r  n  j,  m- 1 

=  [.2  p  .x  ]  .2  Y  i  .x 
i=o  ^n-i  j=o  \m-l-j 


j+1 


m+n-1  min (1, m-1) 

0  2  .2 

l=o  j=o 


P  n  ,  -Y  i 
n-l+ j  rm-l- j 


1+1 


Hence  we  can  write  (3.^.5)  as 


m+n-1  ,  min(l,m-l) 

l(x)  =  „  2  x^  1  [  .2  |3  0A_  .q  .  n 

v  '  l=o  j=o  n-l+ j+n-j-l 


pn-l+ j^.m-1-  j  ^ 


n 

+  .2  (3  .x' 

i=o  n-i 


. 


I'  V }■( 
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At  this  stage  we  find  it  helpful  to  write  out 
a  few  coefficients  of  this  polynomial.  Some  coefficients 
are : 


xra+n  : 

q  P 
Ho  0 

- 

^0^0 

m+n-1 

X 

o 

CO. 

1 — 1 
o' 

+ 

qoPl  - 

pn  y 
^1 r  o 

m+n- 2 
x  : 

q2Po 

+ 

ql^l  + 

qoP2 

Pen 

p2Yo  *  Phi  -  P(h2  ' 


It  is  now  evident  that  these  coefficients  form 
a  triangular  matrix  if  we  regard  the  f3^  and  the  y^  as 
unknowns  and  if  we  equate  these  coefficients  to  any  chosen 
polynomial  provided  that  not  all  the  relevant  coefficients 
p.,  q.  are  zero. 


For  example  *  suppose  that  o  =  m  +  n  +  1.  Then 
let  us  try  to  make  ¥(x)  assume  the  form 


q(x) 


(x-x2 


m+n 

x 


)(x-x3)...(x-xm+n+1 
-  (x2+...+xm+n+1)  X 


) 

m+n- 


1 


We  then  have  that 


qoPo  -  PoU  =  1 

qlP0  *  qoPl  -  Pho  -  poU  =  -(x2+x3+--'+Xm+n+l) 


and  we  can  determine  values  of  f3Q.,  yq  to  satisfy  these 


. 


...  V;  ...  ' 
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equations  unless  pQ  and  qQ  are  both  zero.  Similarly  we 
may,  in  this  example,  make  Y(x)  assume  any  of  the  forms 
¥^(x),  k  =  1,  ...  ,  m+n+1.  This  would  force  us  to  conclude 

that  all  the  A^  are  zero  and  by  Lemma  1,  Q(x)  would  not 
be  the  best  approximation.  Hence  we  are  forced  to 
conclude  that  pQ  =  qQ  =  0  in  this  example. 

We  can  extend  this  argument  to  show  that 
PQ  =  P1  =  P2  =  . . .  =  Pd-1  =0  d=m+n+2-  o 

qo  =  ql  =  q2  =  *  *  *  =  qd-l  =  0 

if  there  are  o  extrema,  which  completes  the  proof  of 
Chebyshev‘s  theorem. 

It  may  be  noted  that  the  proof  of  the  above 
theorem  does  not  require  that  the  extrema  alternate  in 
sign.  It  seems  that  Chebyshev  was  aware  of  this 
property  of  the  rational  function  of  closest  approximation; 
but  it  is  not  explicitly  stated  in  his  longest  paper  on 
this  topic  (see  [2]),  and  we  have  been  unable  to  trace 
an  explicit  statement  of  it  in  the  recently  published 
Ouevres.  Priority  in  the  statement  of  this  property  is 
a  question  of  historical  interest.  It  may  be  noted 
however  that  Achiezer*s  statement  of  the  theorem  complements 


. 

■ 
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the  above  proof  since  it  shows  that  the  extrema  of 
r(x)  =  Q(x)  -  f ( x ) 

must  alternate  in  sign  when  Q(x)  is  the  function  of 

closest  approximation  among  the  set  of  functions 

p  (x) /q  (x)  . 

*n v  ' '  hn v  ’ 


Section  3.5  Rational  Approximations  to  Sets  of  Ordinates 


In  this  section  we  shall  attack  the  problem  of 
extending  the  argument  used  in  the  Chebyshev  concepts  of 
chapter  il  section  6  to  the  rational  approximation  case. 
However,  it  will  be  found  that  the  use  of  determinants 
is  extremely  useful,  if  not  necessary.  The  problem  is 
difficult  and  only  the  preliminary  work  and  a  start  on  the 
solution  along  with  a  few  examples  will  be  found  here. 


We  will  investigate  the  general  approximant 

.n 


(3.5.1) 


y  = 


an  +  3-n-l^  •  •  «  +  &ox 


bn  +  Vlx  +  •••  +  box 


n 


For  definiteness,  we  shall  base  most  of  our  arguments  on 
the  very  simple  approximant 


(3.5.2) 


y  = 


ai 


+  aQx 


b.,  +  b  x 
1  o 


but  most  of  our  arguments  will  apply  with  little 
modification  to  the  general  approximant.  Indeed  we 


shall  explicitly  point  out  any  limitations.  Notice 


i  o 


II 

' 

i  ca  :  ■  .....  10  Dns  ;j  [ t/b  '■  V:: lb 


. 


' 


that  no  proper  rational  approximation  can  be  simpler 
unless  it  reduces  to  a  polynomial. 
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We  can  write  (3.5-2)  as 

(3.5.3)  ax  +  aQx  -  b1y  -  bQxy  =  0 

or  as 


3-5.l0  [1  x  y  xy]  |~  a± 

a 

-b 


Notice  that  (3.5.2)  has  only  three  free  constants 
the  points  (x±,  y±),  1=1,  2,  3,  lie  on  (3.5.2), 
have  from  (3-5.^) 


(3.5.5) 

'  1 

X1 

yl 

xiyi 

a  -i 
-L 

1 

X2 

y2 

x2y2 

a 

0 

1 

x3 

y3 

x3y3 

~bl 

1 

X 

y 

x  y  _ 

-b 

L  O  - 

the  last  equation  of  (3.5-5)  expressing  the  fact 
(3.5.^)  is  true  for  any  (x,  y)  on  the  curve. 


.  If 
then  we 

0  1 
| 

0  i 
0 

0  J 

that 


The  determinant  of  (3-5-5)  must  vanish.  Hence 
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the  equation  of  the  curve  of  form  ( 3 . 5 -2) through 


(x^  y±). 

i  =  1,  2, 

3, 

is 

(3-5.6) 

1 

X1 

yi 

xiyi 

1 

X2 

y2 

x2y2 

=  0 

1 

x3 

y3 

x3y3 

1 

X 

y 

x  y 

Clearly,  our 

reasoning 

extends 

to  higher 

order  approximants . 

Now  our  aim  is  to  approximate  an  arbitrary 
continuous  function  by  a  suitable  rational  function; 
but  this  problem  is  much  too  difficult  to  encompass  in 
a  single  step.  We  begin  with  a  much  simpler  problem. 
To  approximate  2n  +  2  ordinates  by  a  rational  function 
of  the  type 


(3.5.7) 


y  - 


a^  +  a^  ,x  +  . . .  +  ax 
n  n- 1  o 


n 


b  +  b  x  +  .  .  .  +  b  x 
n  n-1  o 


n 


Since  the  degree  of  the  polynomial  equation  for  p 
is  one  greater  than  the  degree  of  the  denominator,  the 
odd  degree  denominator  yields  a  more  difficult  problem 
in  some  respects.  Hence  we  may  as  well  face  at  once 
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the  difficulty  of  an  odd  degree  denominator. 

In  keeping  with  our  purpose  of  simplicity  of 
presentation,  we  shall  base  most  of  our  arguments  on  the 
simple  approximant  (3.5*2),  which  corresponds  to  n  =  1. 

Suppose  we  use  (3-5.5)  to  approximate  four 
ordinates  y^,  i  =  1,  2,  3*  4.  Obviously  we  cannot  fit 
them  exactly  by  a  curve  of  the  form  (3.5-2).  Therefore 
let  us  use  the  notation 

y±  +  pf±(v±),  i  =  1,  2,  3,  4, 

for  the  ordinates  of  the  approximant. 

a-,  +  a  x 
1  o 

y  =  - - 

b-,  +  b  x 
1  o 

at  the  points  x^,  i  =  1,  2,  3)  4.  The  functions  f^(v^) 
are  continuous  differentiable  functions  of  their  arguments 

4- 

and.  vary  continuously  between  -  1.  We  wish  to  emphasize 
by  the  notation  that  the  results  we  shall  obtain  are 
quite  general  and  in  no  way  restricted  to  any  choice  of 
function . 


Since  the  maximum  deviation  between  the  ordinates 
and  the  approximate  ordinates  is  -  p^,  at  least  one  of  the 
fi(vi)  must  be  -  1. 


From  (3.5-6), 
subject  to  the  equation 
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the  ordinates  y^  +  f^(v^)  ane 


(3.5.8) 


F(p ;  v±) 


1 

X1 

yi+ffvfp 

x1(y1+pf1(v1)) 

1 

X2 

y2+f 2( v2)p 

x2(y2+pf2( v2) ) 

1 

x3 

y3+f3^  VP 

x3(y3+pf3( v3) ) 

1 

x4 

y4+f4(v4)p 

*4(y4+Pf4(v4)) 

We  now  seek  to  minimize  p  with  respect  to  the  v^.  The 
necessary  conditions  are  of  the  form 

(3.5.9) 


1 

X1 

pf{(v1) 

x1(y1+pf1( v±) ) 

1 

xi 

yi+pfi( vi) 

px]f1I(v1) 

1 

x2 

0 

x2(y2+pf2( v2) ) 

+ 

1 

X2 

y  2"^p  ^2  ( v2 ) 

0 

1 

x3 

0 

x3(y3+pf3( v3) ) 

1 

x3 

y3+pf3(v3) 

0 

1 

x4 

0 

x4(y4+pf4(v4)) 

l 

x4 

y4+pf4( V4) 

0 

where  we  have  differentiated  with  respect  to  v^.  Three 
similar  equations  can  be  obtained  by  differentiating  with 
respect  to  and  v4 .  At  a  maximum  we  must  have 

=  0 

8v . 


F  dp  +  F,,  dvn- 
P  K  vq  1 


but  since 


0 
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dp  _  Fvj 


that  is j  we  must  have  PVi  =0,  F  ^  0.  Now 

(3.5.9)  is  satisfied,  by  either  of  the  following 

(a)  f *( vx)  =  0  that  is,  at  the  maximum  or 

minimum  value  of  f-^(v^) 

(b)  f^(v^)  ^  0  but  sum  of  determinants 

is  zero; 

and  similarly  for  the  other  three  equations  of  this 
type . 


Let  us  take  hypothesis  (b)  first.  If  we 
multiply  the  equations  of  type  (3.5-9)  by  fq(  vj_)/pf  vj_) 
and  add  the  four  we  obtain 


(3.5.10) 


1 

X1 

fi<v 

^(yypqtvp) 

1 

X1 

yi+pf1( Vq) 

x1f1( v±) 

1 

X2 

^2  ^  V,2  ^ 

x2(y.2+pf2(  v2) ) 

+ 

1 

X2 

y2+pf2( v2) 

x2^  2 ( ^2  ^ 

1 

x3 

f 3<  v3> 

*3(y3+pf3( v3) ) 

1 

x3 

y3+pf3(v3) 

x3f3( v3) 

1 

x4 

X4(y4-tpf4(v4)) 

1 

x4 

y4+pf4( v4) 

x4f 4( v4) 

The 

L.  H.  s.  of  (3.1 

5.10) 

is 

simply  F  . 

Hence 

hypothesis  (b)  does  not  yield  the  necessary  conditions  for 


a  minimum. 


■ 


■ 


■' 

■  ;  -  ii'1  -  •  vi  0  «  (,Q) 


N 


' 


■  ■’  ■  •  VXi;}  ,V. 

V".,' 


*■ 

■ 
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Thus  our  necessary  conditions  for  a  minimum 


must  be  obtained  from  hypothesis  (a); 
fj(v±)  =0  or  fi(vi) 

We  can  now  write  (3.5-8)  in 

(3.5.11) 


1 

X1 

yq+^P 

0) 

+ 

1 — 1 

>> 

1 — 1 

K 

lP) 

1 

X2 

y2+e2P 

x2(y2+e 

2P ) 

1 

x3 

y3+€3p 

x3(y3+e 

3P ) 

1 

x4 

yq=4P 

x4(y4+e 

4P) 

that  is 
=  ±  1 


the  form 


=  0 


where  the  are  -  1  we  don’t  know  which. 


It  might  appear  at  this  point  that  we  aught 
to  demonstrate  the  existence  of  a  minimum.  However  we 
shall  defer  this  for  the  present  since  we  have  no  know¬ 
ledge  of  the  kind  of  difficulty  which  an  existence 
proof  will  encounter  and  must  surmount.  Actually,  there 
is  a  perfectly  good  existence  proof  which  depends  on 
continuity  considerations  but  we  should  like  to  achieve 
a  proof  which  rests  on  simpler  considerations. 

Let  us  now  turn  to  the  question  of  uniqueness. 
Suppose  we  have  a  function  of  the  type  (3-5.2)  which 


■ 
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approximates  the  ordinates  y^,  i  =  1,  2,  3,  by  the 
ordinates  y1  +  p,  y2  -  p,  y^  +  p,  y^  -  p,  that  is 
suppose  the  deviations  at  the  points  (x^,  y^)  are  of 
equal  magnitude  and  of  alternating  sign. 

Suppose  also  that  the  denominator  b^+bQx  does 
not  vanish  in  the  interval  x-^  to  x^  -  we  assume  that 
the  x^  are  arranged  in  ascending  order  of  magnitude. 

Then  we  assert  there  can  he  only  one  function 
of  the  above  type .  For  suppose  there  are  two  different 
solutions  of  (3.5.11)*  say  p(x)/q(x)  and  p(x)/q(x) 
with  maximum  deviation  p,  p  respectively.  Now  p,  p 
must  be  different  for  otherwise  the  rational  functions 
would  be  equal,  too.  For  x^,  i  =  1,  2,  3*  4,  the 
difference 

P(x)  _  p(x)  =  pq  -  pq 

q(x)  q(x)  qq 

assumes  the  values  (-l)1(p  -  p)  ^  0,  thus  being 

subjected  to  at  least  three  changes  of  sign.  The 
numerator  is  of  degree  two,  and  can  account  for  at  most 
two  sign  changes.  Hence  this  leaves  at  least  one  sign 
change  to  the  denominator;  that  is,  one  of  the  approxi¬ 
mations  must  have  a  pole  within  the  range  of  approximation. 
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The  argument  can  easily  be  generalized. 

We  have  no  reason  as  yet  to  favour  the  existence 
of  one  approximant  of  a  given  type  nor  the  existence  of 
another  approximant  of  this  type.  It  is  sufficient  that 
coexistence  is  denied.  Hence.,  if  we  can  demonstrate  the 
existence  of  any  one  approximant,  the  argument  is  complete. 
In  particular,  if  an  approximation  of  equal  and  opposite 
deviation  exists,  we  can  do  no  better. 

To  demonstrate  existence  algebraically,  we 
shall  have  to  examine  closely  determinants  of  the  type 
(3.5.11) . 


At  this  point  it  appears  helpful  to  work  out  a 
few  examples  to  fix  the  above  concepts. 

Section  3. 6. 2  Examples 

Example  1  To  fit  an  approximation  of  the  type 

a]_  +  aQx 


X 

0 

l 

2 

3 

y 

-l 

0 

2 

_ 

1 

to  the  points 


. 

■  •  '■  ;  ■  1  '■  -  on.js 
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We  will  look  for  an  approximation  of  equal  and 
alternating  deviations,  ^  p.  Employing  (3.5.11)  with 
=  (-1)  and  evaluating  for  p,  we  obtain 

p  =  T,/7/4.  This  would  contradict  our  general 
argument  if  both  values  of  p  were  admissible.  Let  us 
first  take  p  = 


The  curve  through  the  points 


X 

0 

1 

2 

y 

-l  -77/4 

•77/4 

2  -77/4 

is 


1  0  -W7/4  0 

1  1  +\T f/4  nT7/4 

1  2  2-aT7/4  4-nT7/2 

lx  y  xy 


=  0 


that  is 

y  =  3  (3-V7)  x  -  3 

4  (4-77)  +  (V7-l)x 


The  crucial  point  is  the  zero  of  the  denominator  at 


x 


(W7) 

(V7-1) 


which  lies  outside  the  interval  [0,3]-  Hence  the 
approximation  is  acceptable. 


. 
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Let  us  investigate  the  value  p  =  -h/7/4.  The 

curve  through  the  points 


X 

0 

1 

2 

y 

-1W7/4 

-V7/4 

2-K/7/4 

is 


1  0  -1W7/4  0 

1  1  -s/7/4  -nT7/4 

1  2  2W7/4  4+nT7/2 

lx  y  xy 


0 


If  we  evaluate  the  determinant  we  find 


y 


3/7  s 


9 


4  +  n/7  -  (l+s/7)x 


The  denominator  has  a  zero  at 

x  =  - L  ~  l.o 

1  +  ^7 

That  is,  the  approximation  becomes  infinite  within  the 
interval  (0,3)  which  is  in  agreement  with  our  uniqueness 
argument,  since  only  one  of  the  approximations  fulfills 
the  stated  conditions. 


In  this  simple  case  we  have  demonstrated  the 


. 


' 


■  •'  ■  1 


- - - ••• 
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existence  of  an  approximation  of  the  prescribed  type 
with  the  necessary  properties.  Prom  the  denial  of 
co-existance,  we  can  do  no  better. 
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We  are  now  fairly  well  alerted  to  the  difficulties 
which  any  proofs  along  this  line  must  surmount.  However, 
let  us  look  at  another  example  which  exhibits  the  difficulty 
in  a  different  way. 


Example  2  To  fit  the  four  points 


X 

0 

1 

2 

3 

y 

-1 

0 

2 

5 

approximately  by 


y 


ax  +  b 
cx  +  d 


Using  (3.5.11)  with  e ^  =  (-l)i-\  we  obtain 

+  V19 


p  =  +  1  x 


~  +2.09,-0.09 


A  plausible  guess  is  that  the  smaller  absolute 
value  of  p  is  the  one  that  we  want,  that  is 


nT19 

4 


P 


+  1 


■  -  ,  o  : 

1  .  . :■  .  "■  ■.  jri:t  ;vj-;:w 
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The  curve  corresponding  to  this  value  of  p  is 


1  0 

1  1 

1  2 

1  x 


-nT19/4 

-I+n/19/4 

3-n/I9/4 

y 


o 

n/19/4  -l 

6-  -s/19/2 

xy 


=  0 


Thus 


1  -(8Vl9  -  19)  +  (llVl9  -  3D 
^  (8  -  VL9)  -  (5  -  V19) 


X 


X 


7  +  \r19 

where  the  pole  is  at  x  =  J - - — - 

is  outside  the  range  [0,3l. 


11.36 


which 


Let  us  now  test  the  other  value  of  p,  where 
\T IQ 

p  =  +1  +  — .  The  curve  corresponding  to  this  value 

of  p  is 


0 

1 

2 

x 


nT19/4 

-(1W19/4) 

3+  s/19/4 


y 


0 


•(1+49/4) 

6+/I9/2 


xy 


=  0 


Thus 


y  = 


1  (8^19  +  19)  -  ( 11  /l9  +  3l)x 

4  (8  +  <sTi9)  -  (5  +  ^19) 


x 


'■■ft  ;!V  sjrtii  A,r 
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where  the  pole  Is  at  x  =  - — -  ~  1.32,  which  Is 
within  the  range  [0,3].  Therefore  this  approximation 
is  not  acceptable.  Let  us  take  one  more  example. 

Example  3  To  fit  the  four  points 


X 

0 

1 

2 

3 

y 

-1 

2 

-2 

5 

by  y  =  -------  using  (3. 5. 11)  with  e  =  (-1)1"1, 

cx  +  d  1 

we  have  p  =  +  "i  “  "Qp  —  +  1.6,  +  3*4.  Taking  the 

least  absolute  value  of  p  first  we  have  the  curve 
given  by: 


1 

0 

3 

2 

- 

n.3 

0 

1 

1 

- 

1 

2 

+  ^  V13 

1 

2 

+  If  ^13 

1 

2 

1 

2 

-  |Vl3 

1 

-  |^13 

1 

X 

y 

xy 

Thus 

=  1  -8  nT13  +  25  +  (7nT13  -  23)x 
4  2  -  sfl3  -  (3  -  ^13)x 

7+^13 


where  the  pole  is  at  x 
within  the  range  [ 0 , 3 ] 


4 


~  2.65  which  is 
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Now  let  us  examine  the  other  value  of 
where  the  curve  Is  given  by: 

l  o  |  +  TfTl3 

1  1  -  -i  -  3 

1  2  j  +  j  Vi3 

lx  y 

Thus 

y  =  1  (8sTl3  +  25)  -  (7sTl3  +  23)x 

^  nTi3  +  2  -  (nT13  +  3)x 

where  the  pole  is  at  x  =  — — -  ~  0.85  which 
lieswithin  the  range  [  0  ^  3  ]  • 

Thus  both  of  our  approximations  in  this  example 
have  poles  within  the  interval  of  approximation.  Hence 
we  have  been  unable  to  find  an  acceptable  rational 
approximation  of  the  required  form  to  the  given  ordinates. 

We  refuse  to  accept  approximation  with  poles, 
because  we  are  trying  to  approximate  -  in  the  end  - 
finite  continuous  functions. 

Section  3.5.3  Interpretation  of  the  Examples 


0 


1 

2 


Jl3 


i  +  \  nTi3 


xy 


=  o 


In  examples  1  and  2  there  was  only  one  admissible 


. 


, 


' 


approximation  which  is  comforting. 


In  example  3  there  was  no  admissible  approximation 
both  approximations  had  poles  within  the  strip  of  points. 
This  example  is  of  considerable  theoretical  importance 
since  it  denies  the  assertion  that  an  admissible  approxi¬ 
mation  of  the  specified  form  to  a  specific  set  of  points 
always  exists. 

We  can  now  see  what  ingredients  must  enter  any 
satisfactory  proof  along  algebraic  lines:  we  must  be 
able  to  find  a  real  root  of  the  equation  for  p*  say 
F(p)  =  0,  which  is  not  a  root  of  the  denominator*  say 
G(p*x)*  for  x^  <  x  <  xn .  One  way  of  doing  this  would 
be  to  show  that  there  exists  real  values  of  x  such  that 
F(p)  =  (ap  +  b)  G(p*x)*  where  x^  <  x  <  xR  and  a^b 
are  constants.  However *  this  does  not  appear  to  be  easy. 

Example  3  implies  one  other  important  point: 
if  the  ordinates  y-^*  y^,  y^,  y^  are  such  that,,  by  adding 
a  suitable  constant* y *  they  can  be  made  to  alternate  in 
sign*  there  is  no  real  root  of  F(p)  =  0  unless  p  is 
sufficiently  large  to  change  at  least  one  alternation 
in  sign.  Hence*  if  we  suspect  the  existence  of  a  real 
p*  we  can  bound  it  by 


' 
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min  |  y±  +  y  |  <  p  <  max  |  y±  +  y  |  . 

A  possible  generalization  of  this  could  be:  if  P(p)  is 
an  even-order  polynomial  in  p  and  if  the  ordinates., 
y  ,  y^,  • • •  »  ^2n+2  are  suc^  that,  by  adding  a  suitable 

constant,  y ,  they  can  be  made  to  alternate  in  sign,  there 
is  no  real  root  of  P(p)  =0  unless  p  is  sufficiently 
large  to  change  at  least  one  alternation  in  sign. 

Our  examples  do  not  imply  any  contradiction  of 
the  existence  theorem  for  rational  approximation.  Rather, 
they  imply  that  certain  sets  of  abscissae  in  the  interval 
of  approximation  may  be  excluded  in  the  search  for  extrema. 
Werner  has  studied  rational  approximations  of  the  form 
(3-5*2)  to  a  continuous  function  f(x).  He  gives  the 
following  theorem  (see  [24]): 

Theorem  3 -5 -3-1  Let  (x^,  f^),  1=0,  1,  2,  3,  be 
given  and  let  x.,-.  >  x..  The  necessary  and  sufficient 

It  _L  1 

condition  for  the  existence  in  [x  ,  x0]  of  a  continuous 

o'  3 

rational  approximation  such  that 

Q(x.)  -  f.  =  (-l)i[Q(x0)  -  fQ] 

is 

sgn ( f Q  -  f2)  =  sgn(f1  -  f^)  . 


••v  . i-  - ‘  :  -  '■  a  1  ;i  j  :  r.qm£  \ 

1  ->  ■■  -V^.7  ...  -  U'K;  i. 
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It  appears  that  this  theorem  can  be  generalized 
to  include  higher-order  rational  approximations ,  although 
this  has  yet  to  be  done. 

We  will  now  consider  two  methods  for  obtaining 
p.  It  will  be  found  that  the  second  method  gives  p  as 
the  eigenvalues  of  a  symmetrical  matrix  with  the 
associated  eigenvectors  giving  the  coefficients  of  the 
denominator  corresponding  to  that  value  of  p. 

Section  3- 5.^  Methods  for  Obtaining  p 

Let  us  write  the  determinant  (3-5.8)  -  in 
generalized  form  -  as 

(3.5.12) 

i  ->  p  n  m  , 

I  1,  X±,  X±S  .  ..  ,  X.  ,  T}±,  X1T]1,  .  ...  x±  T11  | 

1=1,  .  .  . ,  n+m+2 

where  =  y^  +  (-l)1-1p.  It  is  helpful  to  work  with 

matrices  rather  than  determinants  here.  We  shall  use 
the  notation 

n. 

for  the  submatrix  [1,  x,  . . .  ,  x  ]  of 


order  (m+n+2)  x  (n+l). 
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X  for  the  submatrix  [l,  x.  ...  ,  xm]  of 
order  (m+n+2)  x  (irH-1 )  ^ 

Y  for  the  diagonal  .matrix  of  order  (.m+n+2)  x 
(m+n+2)  whose  diagonal  elements  are  and 

D  for  the  diagonal  matrix  of  order  (.m+n+2)  x 
(m+n+2)  whose  diagonal  elements  are  (-1)1-1  . 

The  above  determinant  (3-5-12)  is  then  the 
determinant  of  the  matrix  [X  .  (Y  +  pD)X  ]  and  can  be 
written  as  |  XnJ,  (Y  +  pD)Xm  |  .  Now  we  use  the  equation 

I  X  ,  (Y  +  pD)X  1  =  0 

as  an  equation  of  determine  p.  Clearly  it  is  equivalent 
to  a  polynomial  equation  of  order  m+1.  Hence  it  would 
be  helpful  to  reduce  the  determinant  to  an  equivalent 
determinant  of  order  m+1.  Two  methods  of  reduction  are 
given  below. 

First  Method 

It  can  easily  be  seen  that  the  set  of  m+n+2 
vectors  XnJ,  DX  are  linearly  Independent.  Let  us  find 
the  inverse  of  [Xn  !  DX^J  partitioning  the  matrix  by 
columns ;  the  inverse  will  be  conformably  partitioned  by 
rows  and  can  be  denoted  by 


' 

, 


(  . 
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Hence  we  have 

[-"]  [xn ;  Dxm]  =  [-5-1-1-]  , 

that  is 

aXR  =  I  aDXm  =  0 

PXn  =  0  |3DXm  =  I  . 


Hence  if  we  premultiply  the  matrix 

[  Xn  !  (Y  +  pD)Xm]  by 


we  obtain 


[ 


_ 

0  '  PYXm  +  pi 


ard  the  determinental  equation 


xn,  (Y  +  pD)Xra  I  =  0 


is  equivalent  to 


PYXm  +  pi  |  =  0 


which  is  of  order  m+1 . 
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This  method  is  easy  to  understand  and  simple  to 
use,  however  it  is  likely  to  be  ill-conditioned  for 
large  n  or  m,  since  Vandermonde  determinants  are  ill- 
conditioned  and  the  determinants  of  the  matrices  Xm  and 
Xn  are  Vandermonde  determinants. 


Second  Method 


We  will  find  it  convenient  to  prove  some  lemmas 


first . 


Lemma  3  The  (p,q)th  element  of  the  matrix 


XmT 


A  is  a  diagonal  matrix 


of  order  (m+n+2)  x  (m+n+2)  whose  diagonal  elements  are 
.  m+n+2  p+q-2 

V  -  15  ill  A  xi  • 


m 


Proof :  We  have  X  =  [1  x  x^  . . .  xm3 .  Hence 


X 


T 


m 


1 

x 

/ 

x* 1 


2 


t  * 


m 


and  since 


A 


0 

0  A' 


0  0 


.  0 

.  0 


A 


m+n+2 
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th  T  p-1 

the  p  row  of  X  A  has  elements  equal  to  x.^  X.. 

^  m  ^  1  1 

Now  X  =  [  1  x  x"1  .  .  .  x  ]  and  hence  the  (p,q) 

element  of  X  ^AX  is 


*ip_1Viq"1  +  x2P'1Axq'1 


m+f  2A  x  P+^2 
1=1  A1  1 


Lemma  4  The  elements  of  diagonal  matrix  A  can  be  chosen 
such  that  XmTAXn  =  0. 

Proof ;  The  lemma  is  equivalent  to  the  statement  that 
i  =  1,  2,  ...  ,  m+n+2.,  can  be  found  such  that 

m+n+2 

2  X  x.  =  0,  s  =  0,  1,  m+n. 

i=l  i  l 

Now  it  is  well  known  that  any  (i+l)  vectors  in  $-space 

n  ~j 

are  linearly  dependant.  Hence  the  ($+1)  vectors  x^  -  s 
s  =  ls  2,  I  are  linearly  dependant  if  x.  /  x.. 

We  shall  suppose  that  x^+^  >  xj_*  The  equations  are 


1 — 1 

1 — 1 

1 - - 

rh 

xi  *  xz+i 

A 

i — i 

i 

i — i 

i 

• 

L  xi  . 

.  A^+l . 

and  can  be  written  as 


A  1 


'* 
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1 

x, 


i+1 


i-1 


x 


ri-\ 

V1 


7+1 


-A. 


1 

x. 


x 


i-1 


which  we  can  write  as 


1 

x. 


x 


i+1 


X 


i-1 


X 


.i-1 

i+1 


-a2/Ai 

-a3/a1 


■VA 


1 

x. 


x 


.i-1 


Using  Cramer’s  rule^  we  can  write 


A 

A- 


k 


1  . 


X,  ,  xn 
k-1  1 


i-1 


X 


k+1 


i_l  i-1  i-1 


x,  x 

k-1  1 


k+1 


1  1 


Xo  X- 


X 


i+1 


X 


i-1  i-1 

X 

3 


X 


i-1 

Li+1 


xi+l 


X 


i-1 
i  +  1 


Now  the  above  determinants  are  of  the  form  of  the 
Vandermonde  determinants  which  can  be  written  as 


1 

1  . 

• 

1 

X1 

X2 

x2 

X1 

xn+l 

x2 

n+1 

xn 

1 

x  n 

2 

xn 

n+1 

whose  value  is 

n+1 

n 

j=2 
i=l 
i<  j 

(V- 

+)• 

Hence  we 

i — 1 
r< 

=  (- 

Dk-2 

Z+l 

n 

1=2 
j=l 
j<i 
j ,  i  A 

(x .  -  y 

v  l 

/  Z+l 

/  11 

J  /  i=3 

/  J=2 
/  J<i 

(x. 


Therefore 


Z+l 

/  Z+l 

n 

i=2 

i — 1 

X 

i 

•H 

X 

'  n 
1=2 

(x 

i^k 

i^k 

Z+l 

'Z+l 

(-1) 

n  (x  - 

i=2 

Xx)  / 

n 

1=1 

i^k 

/g_ j_i 

\  =  (x1  -  xk)-1,  k  =  1,  2,  .  .  . ,  Zl 
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Now  since  x. , ,  >  x.  then  if 
i+l  l 

A-,  is  positive  we  have 
A ^  negative, 

A^  positive, 
and  so  on . 


Hence  for  x^+^  >  x^,  sign  =  -sign(Ai) 

This  completes  the  proof  of  Lemma  4. 

If  the  matrix  [Xn,  (Y  +  pD)Xm]  is  premultiplied 
T 

by  X  A,  where  the  A.  are  chosen  in  accordance  with 
J  m  9  i 

Lemma  4,  we  have 


(3.5,13)  xmTA[V  (y  +  P3)xm] 

=  [0,  X  TAYXm  +  pX  TATX  ] 

L  J  m  m  r  m  m J 

and  the  determinantal  equation 

j  X  ,  (Y  -1-  pl)X  I  =  0 

1  n’  v  r  '  m  1 


is  equivalent  to 


(3.5.14) 
Since  (D) 


X  TAYX  1-  pX  tadx 
m  m  r  m  m 


=  0 


.  .  =  (»l)1_16..,  where  6.  .  is  the  Krone cker 

ij  ij  iJ 


delta,  we  have  from  Lemma  4  that  the  diagonal  elements 


of  AD  are  positive.  Hence  the  second  of  the  above  matrices 
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is  symmetrical  and  positive  definite  while  the  first  is 
clearly  symmetrical.  Thus  we  can  write 

T  T 

X  iADX  =  LL1 
m  m 

where  L  is  real.  Premultiplying  (3.5.1^)  by  L-1  and 

T  —  1 

post-multiplying  by  (L  )  ,  we  have 

(3.5.15)  |  L*1XmTAYXm(LT)"1  +  pi  |  =0. 

All  the  eigenvalues  of  the  above  matrix  are  real  since 
it  is  symmetric.  Also.,  its  eigenvectors  are  orthogonal. 


Denoting  the  eigenvectors  of  the  matrix 

corresponding  to  the  determinant  (3-5.15)  by  r  ., 

J 

j  =  1,  ...  ,  m+1,  the  eigenvectors  of  the  matrix  of 
(3.5.1^)  are  given  by 


tT 

D  =  L  qj 

T 

and  the  orthogonality  relation  r^.  r ^  ==  0  is 

T  T 

equivalent  to  LL  q^  =  0  for  j  ^  k.  But, 


(3.5.16) 


m  m 

VLL\ 


m  m 

q,  X  iADX  q, 
^k  m  m 


m+n+2 

^  I  h  I  qk(xl)qj(x1) 


=  0 


where  q.  (x)  =  a  .  +  a  ,  ,  x  +  .  .  .  +  a  ,  x 

^kv  ’  m,k  m-l,k  o,k 


m 


■  ; 


' 

1  O’ 
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Now,  from  (3.5.16)  It  is  impossible  that  q^(x)  and 
q^(x)  should  both  be  one-signed  at  each  of  the  points 
x . ,  i  =  1,  2,  ...  ,  m-l-n+2 ;  for  if  so,  the  L.H.S.  of 

(3.5.16)  would  be  positive  or  negative.  Hence,  only- 
one,  at  most,  of  the  q^  vectors  can  correspond  to  a 
denominator  of  Q(x)  which  is  one- signed  throughout  the 
interval  of  approximation. 
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CHAPTER  IV 


METHODS  OP  FINDING  RATIONAL  INTERPOLATION  FORMULAE 


Section  4.1  Introduction 

Although  this  thesis  Is  mainly  concerned  with 
"best-fit"  -  In  the  Chebyshev  sense  -  rational  approxi¬ 
mations,  this  chapter  is  devoted  to  a  few  of  the 
algorithms  used  to  find  rational  interpolation  formulae, 
since  there  is  still  a  definite  need  for  interpolatory 
rational  approximations.  The  algorithms  will  be  preceded 
by  a  brief  discussion  of  the  relationship  between 
continued  fractions  and  rational  functions .  The  algorithms 
given  below  are  those  that  give  interpolatory  rational 
approximations  in  either  continued-fraction  form,  or 
proper  rational  function  form  -  that  is,  a  rational 
function  with  a  polynomial  denominator. 


Section  4.2  Some  Relationships  between  Finite  Continued 


Fractions  and  Rational  Functions 


The  form  of  a  typical  continued  fraction  is 


(4.2.1) 


y 


a 


1 


+  b 


2 


a 


3 


+  .  .  . 


which  is  commonly  written  as: 


(4.2.2) 


y 


a  +  DU  +  iaj  +  iaJ  + 


Ip  ■  £ 1  o  Ou  y  .■  : 

y  -■  I  i,  t  Si-  ■  '  .  ' 


c £ ..  ■■"■■■■  i 

! 

•  '  ■■■.  r.:-> .? 

■-  Le  :.qvx  ife  'io  n>c: o'i  -eJ!“ 
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A  convenient  method  of  examining  (4.2.1)  or 
(4.2.2)  is  to  express  the  continued  fraction  as  a 
rational  function,  where  the  successive  higher-order 
rational  functions  are  called  convergents.  The  first 
few  convergents  are: 


a  an  +  bn 
o  1  1 


q- 


p2  =  aoala2  +  ao^2  +  a2bl 
ala2  +  p2 

Examining  the  mode  of  formation  of  these 
rational  expressions  we  see  that  successive  convergents 
can  be  obtained  from  the  difference  equations 


(4.2.3) 


n 


q 


n 


a  p  n 
n*n-l 


a  q 
nm-1 


+  b  p 
+  b  q 


n- 


n- 


2 

2 


If  we  replace  the  b^  by  c^x-x^)  we  see  that 
for  odd  n  the  n"^  convergent  has  the  same  power  of  x 
in  the  numerator  as  it  has  in  the  denominator,  and  for 
even  n  the  degree  of  the  numerator  is  one  greater  than 
the  degree  of  the  denominator.  Hence  we  can  write 


)j r'y'lh  :  .  ic  1  ,  ,  "  :  ::T  J  f:-;  A 

-A  t.*v  .  a  :>L  a  -v„:J  :  :  :  Aw  .  .  ./A  r'n  v/  ■  :  s.  j 


(  )  Tjcf  r  A  '/h  ,  ,Vf  TI 


e/A!  n  o  %o A 
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(4.2.4) 


rnW 


an  +  aix  +  ...  +  anxP 

— - A - 2 - pT I  -  "  =  2P 

pQ  +  exx  +  . . .  +  pp_lXp 


a  +  oux  +  .  .  .  +  auxP 

o _ 1 _ p_ 

+  P2X  +  • • •  +  PpxP 


n  =  2p+l 


Since  the  numerator  and  denominator  of  either 
form  of  (4.2.4)  can  be  divided  through  by  any  one  of  the 
non-zero  coef f icients ,  the  first  form  involves  2p 
independent  parameters  and  the  second  form  2p+l  such 
parameters,  so  that  in  either  case  n  independent 
constants  are  available  for  the  determination  of  the 
approximation.  Hence  the  n^h  convergent  gives  an 
approximation  to  f(x)  by  a  rational  function  of  x  which, 
in  general,  agrees  with  f(x)  at  the  n  points  xQ,  x^,  ..., 

xn  ^  if  a  0  and  if  all  preceding  a's  are  finite. 

The  uniqueness  of  a  rational  function  of  a 
given  form  passing  through  the  points  (x^,y^), 
i  =  1,  2,  . . .  ,  n,  is  easily  shown.  We  proceed  thus: 

Let  us  denote  the  rational  function  passing  through 
the  given  points  by 


-  •'  .  :  ",  '  •'  ■'  - . .  •  -•  •  •  -  ■  ;• 


i  (  .  --)  ar  'I.  i  i!,:  i  r.  ;  i  HC  • 

•'  o  iq  p\»  .  r.;','Ojrf3  ,•  £■..  ,  ,  ,,  ,V  ,1  ■■=  r 
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Suppose  that  there  is  another  rational  function: 

Pq(x) 

W 

such  tha  t 

Pn ) 

-  =  Ya,  i  =  1.  2,  ...  ,  n. 

5n(xi) 

Then  ^(x)  =  pn(x±)  qn(x±)  -  Pn(xj_)  qn(x±)  =  0, 

i  =  1,  2,  .  .  .  ,  n.  (We  are  assuming  here  that 
qn(Xi)  ^  0,  q  (x^)  ^  0.)  That  is,  the  polynomial 
if/  (x)  vanishes  for  n  values  of  x.  Whether  n  is  even  or 
odd,  the  degree  of  if/{x)  is  n-1.  Hence  if/(x)  must  be 
identically  zero.  Therefore  pn(x)/qn(x)  and 
Pn(x)/qn(x)  are  identical  except  perhaps  for  a  factor; 
that  is  pn(x)/qn(x)  may  be  of  the  form: 

A ( x ) Pn ( x ) 

A(x)qn(x) 

However,  we  may  find  that  the  rational 
function  -  of  a  specified  form  -  which  is  to  pass 
through  the  n  points,  has  a  pole  within  the  interval  of 
interest,  or  leads  to  an  inconsistent  set  of  equations 


.  X  ttO  z:-§.i  1  ■'  Itol  ri 3  JTf£V  ( X  ] 

'X  '  :■  ;,jv\  •'  q.  r  . 

:T  1  -2i  (3  i  (IB  ■  ■,  o,  '  3.1.  "L 


3  '.3(  ?J  &  1 
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for  the  determination  of  the  coefficients  of  the  function. 
For  example  finding  the  rational  function  that  passes 
through  (-1,-1)  and  (1,1 )  and  has  the  form 

a 

x  +  b 

we  obtain  a  =  1,  b  =  0  and  the  function  has  a  pole  at 
x  =  0.  If  we  take  the  function  to  be  of  the  form 

a 

bx  +  1 

we  obtain  the  equations 

a  -  b  =  -1 

a  -  b  =  1 

which  are  inconsistent. 

Section  4.3  The  Inverse  Divided  -  Difference  Interpolation 

Formula 

We  can  consider  Newton's  divided-difference 
polynomial  interpolation  formula  (Chapter  2,  section  4), 
with  an  error  term,  as  the  identity  which  results  from 
writing 

=  uQ(x) 


(4.3.1) 


f(x) 


'  ■:  "  ■  t'  l-  .  ,  1  /  T  >  -fc.  ,f:a£nj.:j'tn  j.:;  "  b  t-ii.r  ool 

■ 

'■  •  i  .  :  : 

■  ike'  '  ■b-  ....  ..  :  ■  ....  1  '  :  O',  r.  . :  OC:  O.  ■■■■•  ■„  elv 

1  "  •  ■  '  t  b  nv::  •.*:  ,  <  Oy 
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and  performing  the  successive  substitutions 

(4.3.2)  uk(x)  =  uk(xk)  +  (x-xk)uk+1(x), 

k  =  0,  1,  .  .  . ,  n-1. 


with  the  abbreviation 

(4.3.3)  uk(x)  =  ftx0^  x!^  •••  *  xk-l3  ' 


The  algorithm  for  the  calculation  of  the 
successive  divided  differences  follows  directly  from 
(4.2.2)  and  (4.2.3).,  with  x  =  xk^  in  the  form 

(4.3.4) 

f[V  x1, . . xk_2,  xk_1>  xk]  = 


f  [x  X-,  ,  ...  ,  X,  n  J  X,  ]  -  f  [x  J  ,  X,  o  J  x.  ] 

L  O'*  1*  3  k-23  kJ  L  o3  3  k-23  k-lJ 


x,  -  x,  n 
k  k-1 


The  result  of  assuming  that  the  (n+l)^^  divided  difference 

t  h 

un+i(x)  is  identically  zero  (or  that  the  n  divided 
difference  is  constant)  is  the  equation  of  the  polynomial 
y(x),  of  degree  <  n,  which  agrees  with  f(x)  at  the  n+1 
points  xq,  x^  ...  ,  xn .  If  un+1(x)  actually  vanishes 
identically,  then  y(x)  is  identically  equal  to  f(x). 


in'  ,  .  : 
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A  variety  of  other  identities  may  be  obtained 
in  a  similar  way  by  making  use  of  other  sets  of  trans¬ 
formations.  In  particular,  the  substitution  sequence 


(4.3.5) 


f(x)  =  vo(x) 


vk(x) 


w  + 


x-xk 

vk+iM 


k  =  0,  1,  2,  .  .  . 


leads  to  the  following  result 


(4.3.6)  f(x) 


x  -  x. 


V  (x  )  +  - - — 

0  °  I  vl(xx) 


+ 


x  -  x1  j  x  -  x§2 
|  v2(x2)  +  I  v3(x3) 


if  truncated  after  three  terms.  Thus,  more  generally, 
we  are  led  to  the  continued-fraction  representation 


(4.3.7)  f(x) 


a  + 
o 


X-X, 


X-Xj  1 

I  a2 


+ 


x-x2  I 


+  . . . 


2 


where  a,  =  v.  (x,  ),  and  where,  when  the  fraction  is 
terminated  after  n  divisions,  the  constant  aR  is  to  be 
replaced  by  aR  +  (x-xn)/vn+1 (x)  in  the  last  denominator. 
If  we  then  set  x  =  x^.,  0  <  k  <  n,  the  fraction  terminates 
before  the  residual  (x-xn)/vn+^(x)  is  introduced.  Since 
(4.3-7)  is  an  identity,  the  result  of  replacing  l/vn+1(x) 
by  zero  -  that  is,  terminating  the  fraction  with  aR 
will  give  a  function  Rn+1(x)  which  agrees  with  f(x)  at 


/;■  .1  '  '  '  -  1  .  ."r,  ;■"!  /  / 

,  j  5  0  f  S  :  I  t'S  E?  £ 

■■  ■  :  '  ■  ' 


eiaffl  :  ; r:' v, .'t  “l 


'.  1  ..? 

,  J.1  > 

.  o.;'  1..  xnJ  :,6-  ,  gr  . j- .■  ? 

,  "  y  xh  /  (..  '  .  f  i.  '‘L1  Xl:i';  ■'  v.  1 1  i: 
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the  n+1  points,  xq,  x^,  ...  ,  xn,  under  the  assumption 

that  the  constants  a  ,  a^  .  ..,  an  actually  exist  and 
that  the  portion  of  the  truncated  fraction  preceding 
x-x^  does  not  vanish  when  x  =  x^,  k  =  0,  1,  ...  ,  n-1. 

If  we  introduce  the  notation 


(4.3.8) 


a. 


k 


reference  to  (4.3.5)  gives 


zq[x]  =  f(x),  z1[xo,x] 


z2(x0,xrx] 


z1[x0,x]-z1[xQ,x1] 


and  in  general 


(4.3.9) 


Accordingly,  we  have 


(4.3.10)  Z]^ [ Xq ,  x1 


xk  -  xk-l 


X,  X,  ]  -  z,  i  [  X  , 

k-2*  kJ  k-lL  o" 


•  •  y 


xk-2,xk-l-* 


4 
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Thus  the  z,  [x  ,  .  .  . ,  x,  , ,  x,  ]  is  the  kth 

kL  o*  k-1  k 

inverse  divided  difference  of  f(x)  relative  to  xk  ^ 

and  x^.  We  shall  refer  to  the  quantity  defined  by  (4.3.10) 
t  hi 

as  a  k  inverted  difference  of  f (x) .  The  inverted 

difference  (4.3.10)  is  symmetrical  in  its  last  two 

arguments  x^  ^  and  x^.  It  is  not  generally  symmetrical 

in  its  other  arguments .  Hence  it  must  be  formed  from  the 

specific  inverted  differences  z,  [x  ,  .  .  .  .  x,  0.  x,  ,  ] 

^  k-1  o  k-2*  k-1 

and  ^[xQ,  ...^  x^  2,  x^]  which  possess  their  first 

k-1  arguments  in  common.  Hence  the  following  calculational 
arrangement  is  most  convenient: 


X 

f(x  ) 

o 

1 — 1 
X 

f  (Xf) 

Z1[VX1] 

x2 

f  (x2) 

z1[x0,x2]  z2[xo,x1,x2] 

x3 

f(x3) 

z1[x0,x3]  z2[xo,x1,x3] 

where  the  diagonal  elements  are  the  desired  constants 
aQ,  a^,  a.^,  ...  which  appear  in  (4.3.7) 

Section  4.4  Thiele^  Continued-Fraction  Expansions 

t  h 

Whereas  the  k  inverted  difference., 

z,  [x  ,  .  .  .  ,  x,  x,  ,  j  x,  ]  ,  of  a  function  f  (x)  was 

ku  o’  k-2  k-1*  k  v  ' 


' 

•  -r/--  J  r ... 

1  '  B,  L  {  .  ,  4:,  ,  ,  . . ;  r. 
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symmetric  in  only  its  last  two  arguments,  it  happens 
that  the  quantity 


(4.4.1) 


PkCV-—  xk] 


[k/2] 

2 

$=o 


Jk-2i 


[x. 


'  xk-2iJ 


is  symmetrical  in  all  of  its  k+1  arguments.  Here  the  last 
term  on  the  right  is  z0[xQ]  if  k  is  even,  and  is 
z^[x  ,x^]  if  k  is  odd.  This  quantity  is  often  known  as 
the  kth  reciprocal  difference  of  f (x)  . 


In  particular,  we  have 


p  [x  ]  =  z  [x  ]  =  f(x  ) 

*o  o  o  o  v  o ' 


Pl[Xo'Xl]  =  ZdXo<Xl]  = 


xX  -  xo 


1L  O'  1‘ 


f(xx)  -  f(x  ) 


and  calculation  shows  that 


P2fxo’xl'x2]  = 


z  [x  ]  +  z0[x  ,xn ,xn] 
oL  oJ  2L  o'  1'  2 


x  f  (fn-fn)  +  xnf-,(f0-f  )  +  xnfn(f  -  f ) 
o  ov  1  2'  1  2  o'  2  2 v  o  ±J 


xo(frf2)  +  x^Vh)  +  x2<fo-fl) 


in  which  case  the  symmetry  is  apparent. 

Since  (4.4.1)  implies  that 

Pk[xo---'xk]  -  Pk-2[xo---’xk-2] 
=  zk[xQ, ...,xk], 


(4.4.2) 
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reference  to  (4.3.9)  shows  that  the  successive  reciprocal 
differences  may  be  obtained  by  use  of  the  recurrence 
relation 


(^. 4. 3) 

Pk[v-’xk]  = 


X,  -  X, 

k  k-1 


pk-l*-xo'  ’  ’  *'xk-2'xk]  "  pk-l^xo* *  *  * jXk-2jXk-l^ 


+ 


pk-2 


•,Xk-2 


] 


While  this  formula  is  less  simply  applied  than  (4.3.9)* 

t  h 

the  symmetry  of  the  k  reciprocal  differences  permits 

t  h. 

its  calculation  from  any  two  of  the  (k-l)  reciprocal 
differences  having  k-l  of  its  arguments  in  common,, 
together  with  the  (k-2)  reciprocal  difference  formed 
with  those  arguments.  Hence  a  reciprocal-difference 
table  may  be  constructed  in  the  convenient  form 


X 

f(x  ) 

0 

p! (x0^x1 ] 

X1 

f  (xx) 

p2[xo,x1,x2] 

X2 

f  (x2) 

P1[x1,x2] 

P2[x1,x2,x3] 

x3 

f(x3) 

p1[x2,x3] 

Prom  this  table  we  may  determine  the  coefficients  of 
(4.3-7)  by  combining  (4.3.8)  and  (4.4.2)  -  a  procedure 


'■  VV .  ''  \,h.  j,  :'V',  :  ; 

■  '  '  -V  X  .0 
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due  to  Thiele  -  so  that 


a 


o 


f(x0)'  al  =  Mxo,xl^a2  =  MXo'xl'x2^"f (xo) 


P3[xo,x1,x2,xc,]  -  p1[xo,x1] 


and  so  forth.  Hence  the  required  coefficients  are  formed 
from  reciprocal  differences  appearing  in  the  forward 
diagonal  beginning  with  f  (x  ) .  Because  of  the  symmetry, 
the  data  from  the  same  table  are  available  for  the 
determination  of  formulae  in  which  the  ordinates  are 
introduced  in  other  arrangements,  by  choosing  different 
paths  made  up  of  suitable  contiguous  diagonal  segments . 

Each  such  expansion  is  identical  with  the  one  that  would 
be  obtained  by  the  use  of  the  inverted-difference  array 
corresponding  to  an  appropriate  reordering  of  the  abscissae, 
but  only  one  array  of  reciprocal  differences  is  needed 
for  the  formation  of  the  entire  set.  Thus  the  use  of 
reciprocal  differences,  rather  than  inverse  differences, 
is  generally  advantageous  only  if  several  such  formulae 
are  required. 

Section  4.5  The  Error  Term  of  Thiele  Expansions 


Let  us  examine  the  difference  between 
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(4.5.1) 


y(xQ,x1, . . .,xn,x) 


x-x. 


y  + 
J  o 


+ 


x-xi  | 

I  a2 


+  .  .  , 


x~xn-2  1 

I  an-l 


x-x 


+ 


n-1 


n 


which  passes  through  the  points  (x  ,y  ),  i  -  0,  1,  2,  .. 
n,  and 

(4.5.2)  y(xo,x1, . . .,xn_1,x,x)  = 


x-x. 


x-x- 


y  + 

J  o 


a. 


+  m  +  X~Xn-2  1  +  x-xn-l 


ln- 1 


an(x) 


where  we  use  the  notation  a^(x)  to  remind  us  that  (x  ,  y 

n '  '  v  n  J  n 

must  be  replaced  by  (x,y) . 


Now.,  equation  (4.5.2)  is  of  no  practical 
value  for  it  shows  merely  that  we  can  terminate  the 
continued  fraction  to  obtain  an  identity. 

To  obtain  an  error  term  we  must  find  an  expres 

sion  for 


pn+]>> 

Vn(x) 


pnw 

qnU) 


From  difference  equations  (4.2.3)  we  have 
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(4.5.3) 


P 


n+l 

%+l 


Ml 


a  Ln  p  +  b  .  n  p  ,  p^ 

n+l*n  n+l*n-l  *n 


an+l%.  +  ^n+l^n-l 


q 


n 


b 


n+l 


P  n  q  -  p  q  -i 
*n-l+i  mimi-I 

qn+iqn 


=  -b 


q 


n 


-1  ,pn  Pn-1 


(+  - 


n+1  ln+1  %  qn-l 


Now,  since 


£o 

<J0 


0 

1 


p- 


q- 


a  an  +  bn 
o  ]  1 


we  have 


Pi 

^1 


_2 

qQ 


9l 


Using  (4.5.3)  we  have  in  order 


(4.5.4) 


Pr 


q2 


Pi 

qi 


qxq2  1  ~ 


Pr 


q 


3 


P  r 


q,- 


q2q3 


b-j  b^b^ 


n+l 


Pn 

^n 


=  (-D 


n 


b,b0  .  .  .b  . 

1  2  n+l 


qn+i 


9n9-n+ 1 
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If  b^  in  these  expressions  is  replaced  by  x-x^  ^  we 
see  at  once  that 


pn+l 

%i+l 


P 


n 


% 


vanishes  at  the  n+1  points  x^  i  =  0,  f  n  -  as 

is  expected. 


Let  Pn/qn  in  (4.5.^)  toe  constructed  from  the 
set  of  points  (x^y^).,  i  =  0,  1,  .  .  . ,  n  and  let 

Pn+l/^n+1  constructed  inom  the  set  of  points  (x^, y^)  , 

1=0,  1,  n  and  (X,f (X) ) .  Then 


(4.5.5)  Pn+1(X)  p„(X)  _  f  _  Pn(X) 


<WX)  qn(X> 


qn(X) 


where 


(_1}n  (X-x0)...(X-xn) 

qn(x)qn+1(x) 


qn+1(x)  =  a^q^X)  +  (X-xJq^UX) 


n+lni 


n'  +1-1 


and  a . ,  =  a  fx  ,x  ,x,,  .  .  ,,x  n  ,X) 

nv  n'  o'  1'  '  n-1'  ’ 


n+1 


X  -  x. 


n 


a(xx,x2,  .  .  .,xn_1,xn,X)  -  a(xn,xv...,xj 


o'  T 


n1 


(4.5.6) 


. 

■  :  ' '  .f.’  '..:-yZ 
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ji  -i  7Xv.  t  ./■  :  e-'Tj'icq  :k/  1  sr- 


u: 


in  (4.5.6)  requires  that 


Since  the  calculation  of  qR+1 
we  know  f(X),  (4.5.5)  cannot  be  employed  in  an  error 
estimate.  In  practice,  we  would  replace  X  by  the  nearest 
tabular  point  in  calculating  (4.5.6)  and  thus  obtain  an 
appropriate  error  estimate  from  (4.5*5). 


Section  4.6  Thacher  and  Tukey's  Algorithm  for  Obtaining 


Rational  Interpolates 


This  algorithm  gives,  in  the  simple  case,  the 
same  rational  function  as  that  given  by  the  inverse 
divided-difference  algorithm  of  section  4.2.  However, 
Thacher  and  Tukey's  algorithm  does  not  require  the  genera¬ 
tion  of  a  large  inverted  divided-difference  table  to  find 
the  a^  of  (4.3.7) ,  and  herein  lies  one  advantage  of  this 
algorithm.  A  more  significant  point  is  that  the  algorithm 
can  be  generalized  to  obtain  rational  interpolants  in 
which  the  degree  of  the  numerator  may  exceed  the  degree 
of  the  denominator  by  k  -  where  k  is  an  integer  -  or  vice 
versa.  The  generalization  may  be  achieved  in  the  follow¬ 
ing  ways : 

(1)  by  using  for  the  b.  in  (4.2.1),  functions 
that  vanish  at  x  =  xj_q 

(2)  by  choosing  appropriate  starting  functions. 


The  algorithm  maybe  derived  as  follows:  The 

convergents,  Pj_ (x)/q..  (x) ,  of  the  continued  fraction 

/  v  x-x0  |  x-xi  I 

f  x  =  a  +T - —*  1 2  +  , - — 1  +  ... 

°  |  a aP 


obey  the  relation  (4.2.3)  with  bq  replaced  by  (x-x^  .); 


that  is 


-H  i  :  '  •  •  •  ■: .  ■  o  (  '  c  *  i:  ’■  :  a )  r:.vnl  -- 


■ 
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P± (x)  =  a1p1_1(x)  +  (x-x1_1)p1_2(x) 


q±(x)  =  a±q1_1(x)  +  (x-x1_1) q1_2 (x)  . 

Hence  if  we  develop  functions  t^(x)  that: 

(i)  are  zero  for  x  =  x^, 

(ii)  obey  the  above  relations,  and 

(iii)  t  ^(x)  and  tQ(x)  are  given  appropriate 
initial  values, 

we  can  find  a^,  i  =  1,  2,  ...  through  use  of 


If  we  give  a^  the  value 


(4.6.2) 


t^(x)  may  be  found  from 


(4.6.3)  t.(x)  =  ait1_1(x  )  +  (x-x1_1)t±_2(x  ) 


and  we  can  show  that  t^(x^)  =  0  as  follows:  Using 

(4.6.1)  and  (4.6.2)  we  have 


+  (x1-x1.1)t1_2(x) 


0  . 


IQ\!  o\c 
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To  obtain  the  a^,  i  =  1,  2,  .  .  . ,  of  the 

inverse  divided-difference  algorithm.,  we  set 

*1  -  XQ  =  (*1  -  xo)(-l) 

1  f(xx)  -  f(x0)  f(xx)  -  f(x0) 

Hence  we  have 

t_1(x)  =  -1  and  tQ(x)  =  f(x)  -  f(xQ)  . 

Employing  (4.6.2)  and  (4.6.3)  we  can  obtain  all  the 
remaining  a^  by  recursion. 

For  the  generalization  of  the  algorithm  we 
shall  follow  Thacher  and  Tukey.  The  series  of  functions 
t^(x)  may  also  be  defined  by 

t±(x)  =  q±(x)  f(x)  -  p1(x)  . 

Provided  that  f(x)  is  finite  and  that  q.(x)  does  not 

J 

vanish^  a  necessary  and  sufficient  condition  that 

(4.6.4)  f.(x)  =  P  j  ( x )  /q  j  ( x )  =  f(x) 

is  that 

(4.6.5)  tj(x)  =  q j ( x ) f ( x )  -  p^(x)  =  0. 

Now  suppose  that  for  some  integers  j  >  0  and  k  >  0, 
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fcj-2(x) 

and  tj  ^(x) 

satisfy  the 

conditions 

(4.6.6) 

=  0 

-k  <  i  <  j-2 

(4.6.7) 

=  0 

-k  <  i  <  j-1 

Then  the 

functions 

(4.6.8) 

Vx)  = 

ajtj-i(x) 

and 

(4.6.9) 

Vx)  = 

tJ-i(x)  + 

x-x,-  1 

J  t.  2(x) 

a  J  ”  d 

vanish  for  x  =  x.  (-k  <  i  <  j-2)  since  both  t.  p(x) 
and  t.  -^(x)  vanish  at  those  points.  (The  simple  case 
of  the  algorithm  corresponds  to  k  =  0  and  the  use  of 
the  above  initial  values  of  t  n (x)  and  t  (x)  to  obtain 
the  same  a^  as  those  found  by  use  of  the  inverse 


divided-difference  algorithm.)  They  vanish  for  x  =  x 
since  t  .  ^(x)  and  (x-x  .  ^  vanish  there ,  and  they  can 


j-1’ 


be  made  to  vanish  at  x  =  x  if  we  give  a  .  the  value 

J  &  J 


(4.6.10)  a.  =  -(x.-x.  ., )  t  .  n(x.)/t.  (x  .) 

v  '  J  v  J  J-1'  J-2 v  j;/  j-lv  j' 


A  possibility  of  trouble  comes  when  q.(x)  has  a  root 
near  one  of  the  base  points.  Under  this  circumstance,, 
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subtraction  of  nearly  equal  quantities  can  lead  to 

appreciable  losses  of  significant  figures.  We  can 

avoid  this  difficulty  by  changing  the  recursion  from 

(4.6.8)  to  (4.6.9)  for  one  or  more  cycles.  An  additional 

possibility  of  difficulty  occurs  if  t ^  ]_(xj)  =  0, 

leading  to  an  infinite  a ^ .  This  can  occur  only  if 

P j _ ]_ ( x j ) /q j _  1  ( x j )  =  f  (xj-)  •  Unless  f  (x)  =  PJ._1(x)/qj._1(x) 

so  that  additional  data  are  superfluous,  postponing  the 

introduction  of  x.  and  f (x  . )  will  resolve  the  difficulty. 

J  <3 

Thus  we  have.,  subject  to  the  conditions  that  no 

q.(x.)  vanishes  for  i  <  j,  an  iterative  method  by  which,, 

J 

given  a  set  of  x^  and  corresponding  finite  f(x^),  and 

suitable  starting  functions.,  t  ^(x)  and  tQ(x),  we  may 

construct  a  sequence  of  functions  t  .  (x)  which  vanish 

J 

for  successively  larger  sets  of  x^ . 


Because  of  the  linearity  of  (4.6.5)  .>  of  the 


definition  of  t.(x),  and  of  the  alternative  forms.,  (4.6.8) 

J 

and  (4.6.9)  of  the  iterative  algorithm.,  the  functions 
q.(x)  and  p.(x)  obey  recursive  relations  of  the  same 

J  J 

form,  and  with  the  same  a  namely  (4.2.3)  with 
bn  =  x-xin_1J  or  if  (4.6.9)  is  in  use,  rather  that  (4.6.8) 


n-1 


(4.6.11) 


x-x 


q  (x)  =  q  (x)  +  - q  (x) 

J  J-i  a j  J-h 
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(4.6.12) 


Pj(x) 


Pj.l(x)  + 


x-x 


LI 


Pj.2(x) 


Since  on  each  application  of  the  Iteration  the 
degree  of  q^(x)  is  no  greater  than  the  greater  of  (l) 
the  degree  of  q  .  ^(x)  and  (2)  one  more  than  the  degree 
of  q  ^  ^(x),  an(^  since  a  similar  relation  holds  for  the 
degree  of  p^.(x),  p  ^(x)  an<^  Pj  t^rie  degree  of  either 

the  numerator  or  the  denominator  of  the  rational  inter¬ 
polate  will  increase  by  at  most  one  unit  for  every  two 
iterations . 


It  is  clear  that  the  argument  would  apply 
equally  well  if  any  other  function  which  vanishes  for 
a  =  a^.  2  were  substituted  for  the  factor  (x-x ^  in 
(4.6.8)  or  (4.6.9).  However,,  except  under  rather  special 
circumstances,  such  as  where  the  function  is  known  to  be 

2  p 

symmetric,  so  that  the  factor  (x  -xj_q  )  would  be  appro¬ 
priate,  such  a  modification  does  not  appear  to  have 
particular  merit  for  interpolation  of  functions  of  a 
single  variable  by  ratios  of  polynomials. 


Another  possibility  for  increasing  the  flexibility 
of  the  algorithm  lies  in  the  choice  of  the  starting  func¬ 
tions.  The  only  necessary  restriction  on  the  starting 
function  is  that  t  (x  )  vanish.  If  tQ(x0)  depends  upon 


. 


'  :  -  ;■  .  'i X  I  •,  . 

) 
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more  than  one  parameter,  we  may  use  the  additional 
flexibility  to  make  t  _(x)  and  tQ(x)  both  vanish  at  k 
additional  points,  x  =  x  x  =  x_2,  ...  ,  x  =  x_^, 
where  the  value  of  k  is  determined  by  the  particular 
starting  functions  selected. 

Ordinarily,  the  full  degree  of  possible 
generality  is  not  necessary.  However,  the  cases  which 
follow  are  important: 

(a)  q-]L(x)  =  qQ(x)  =  1, 

while  p  ^(x)  is  a  polynomial  of  degree  k-1,  which  agrees 
with  f(x)  at  x  =  x.  (-k  <  i  <  -1)  and  pQ(x)  one  of  degree 
k,  which  agrees  with  f(x)  at  x  =  x.  (-k  <  i  <  0) 
and 

(b)  p_1(x)  =  pQ(x)  =  1, 

while  q  ^(x)  is  a  polynomial  of  degree  k-1,  which  agrees 
with  l/f(x)  at  x  =  x.  (-k  <  i  <  -l)  and  qQ(x)  one  of 
degree  k  which  agrees  with  l/f(x)  at  x  =  x.  (-k  <  i  <  0), 
respectively.  Case  (a)  ordinarily  leads  to  a  sequence 
of  ratios  of  polynomials  in  which  the  degree  of  the 
numerator  exceeds  that  of  the  denominator  by  k  and  k-1, 
alternately,  while  case  (b)  leads  to  a  sequence  in  which 
the  degree  of  the  denominator  exceeds  that  of  the 
numerator  by  k  and  k-1,  alternately. 


, ) , 


■ 
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All  the  rational  functions  so  obtained  will 
be  minimal  in  that  (neglecting  the  triviality  of  a 
numerical  multiplier  common  to  numerator  and  denominator) 
there  are  exactly  as  many  free  coefficients  to  be 
determined  as  there  are  points  to  be  fitted.  By  a 
suitable  choice  of  k  and  either  case  (a)  or  case  (b), 
one  can  obtain  a  rational  interpolant  with  any  minimal 
combination  of  degrees  for  numerator  and  denominator. 
Since  all  such  are  unique  (up  to  a  common  constant)  all 
minimal  rational  interpolations  can  be  obtained  via  the 
generalized  algorithm. 

The  starting  functions  for  cases  (a)  and  (b) 
are  clearly  the  interpolation  polynomials  for  f(x)  and 
l/f(x),  respectively,  based  on  the  points  x  .  ..,  x  ^ 
for  p_^(x)  or  q  ^(x),  and  on  the  point  xq  as  well  for 
Po(x)  or  qQ(x) . 

If  the  explicit  form  of  the  interpolating 
function  is  desired,  the  starting  functions  must  also  be 
generated  explicitly.  The  following  algorithm  allows 
this  to  be  done  in  an  iterative  fashion. 

Let  ft(x)  or  i p~{x)  be  the  desired  polynomial 
J  J 

of  degree  j+k,  -k  <  j  <  0,  depending  upon  whether  we  are 


/  '  :j  .  :  . 
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developing  p^(x)  or  Qj(x)  •  Thus  we  must  satisfy  one  of 
the  following  sets  of  conditions 


(4.6.13)  ^j(xi)  =  f(xi)  or 

(4.6.14)  ^j(x±)  =  1/f(x1) 

where  -k  <  i  <  j  <  0 ,  Let  Sj(x)  be  a  polynomial  of 
degree  j  +  k  +  1  which  vanishes  for  the  x^  of  (4.6.13) 

zb 

or  (4.6.14).  Then,  if  ^  ,  (x.)  satisfies  (4.6.13) 

J  -L  1 

or  (4.6.14),  and  s  .  ^(xj_)  vanishes  for  -k  <  i  <  j-1,  the 
functions 


(4.6.15)  Sj(x)  =  (x-xj)sj-_1(x) 

and 


(4.6 . l6) 

±  +  [f  (x  .)  ]  1  -  V'i  1±(x1) 

t*  (x)  =  ip  (x)  +  - 4 - J - s  (x) 

J  J_1  s  .  n  (x .)  J_1 


will  satisfy  the  imposed  conditions  for  x  =  x.  (-k  <  i<  j) . 
If  we  use  as  initial  functions 


(4.6.17)  sk(x)  "  ^x_x-k^ 

and 


+ 


[f(x.k)]±:L 


(4.6.18) 
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we  can  generate  our  Interpolation  polynomial  recursively. 
This  algorithm.  It  may  be  noted.  Is  subject  to  round-off 
error,  and  therefore  Is  not  as  precise  as  some  of  the 
classical  methods  of  generating  Interpolation  polynomials 
explicitly.  It  has,  however,  the  advantage  of  simplicity, 
and  of  adaptability  to  successive  introduction  of  data 
points . 


Once  the  starting  functions,  q  -^(x),  P 
t  ^(x),  q  (x),  pQ(x)  and  tQ(x)  have  been  generated,  the 
generalized  algorithm  proceeds  exactly  as  in  the  simple 
case  given  below. 

A  streamlined  version  of  the  simple  algorithm 
(k  =  0)  can  be  written  as: 

Starting  Values 

-lx  i'  ov  i'  v  i'  v  o' 

for  i  =  1,  2,  . . .  n 

q_x  (x)  =  0,  qQ(x)  =  1 

P_l(x)  =  1,  PQ(x)  =  f(xQ) 

Iteration  (  j  =  1 ,  2 ,  . . .  ,  n ) 


a.  =  -(x.-x.  ,)t.  n(x.)/t.  n  (x  .) 

J  J  J-l'  J-2 v  jJ/  j-lv  j' 


XO  Vi.*} 


1 


io6. 


4)(xi>  =  a/j-l(xl>  +  (xi~xj-l^j-2(xl) 
or*  =  pfj.^Xi)  +  (xi_xj.1)^j_2^xi)/aj 

for  i  =  j+1,  j+2,  ....  n 

where  ^J.(x1)  =  tj(x±),  q.(x),  or  Pj(x) 

*  Use  one  or  the  other  for  t  .,  q.,  and  p.  with  a 

j  j  j 

particular  value  of  j . 


Order  of  Calculation 

a.,  t.(x.)j  q.(x)**,  p.(x)** 

j  iv  iv  '  y  ' 

**  These  two  functions  are  calculated  only  if  the 
rational  function  p(x)/q(x)  is  required. 

The  calculation  of  the  a.  -  using  the  simple 

J 

algorithm  -  is  demonstrated  by  the  following  algorithm. 


1 . 

i  <— 

-1 

2. 

tl( 

-  1 

3. 

t2( 

xi)^- 

f (x±)-f (xQ) 

4. 

i  <r- 

-  i  hi 

5. 

i :  n 

< 

>2. 

6 . 

j  <“ 

— 1 

7. 

t2( 

Xj):° 

- = - >  19 
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8. 

aj 

(xj- 

1-Xj)tl(xJ)/t2(xJ.) 

9. 

i — 1 

+ 

*T— D 

I 

•H 

10. 

1 :  n 

> 

>  16 

11. 

t3  <- 

ajt2(x±)  +  (xi-xJ 

12. 

tl(x± 

)<— 

•H 

x 

OJ 

-p 

1 

13. 

t2(xj_ 

)< — 

-t3 

14. 

i+1 

15. 

Go  to 

10 

16. 

j:n 

> 

— *  21 

17. 

j  < 

j+1 

18. 

Go  to 

7 

19. 

m  < - 

■j-1 

20. 

Go  to 

22 

21. 

m  < - 

n 

22. 

Print 

:  a 

.  for  j< - 1  ( l)m 

J 

23. 

Stop 

24. 

End 

]_)  tl(xi) 
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CHAPTER  V 

ALGORITHMS  FOR  OBTAINING  RATIONAL  APPROXIMATIONS 


Section  5.1  Introduction 


This  chapter  contains  a  few  of  the  algorithms 
for  obtaining  rational  approximations  with  emphasis  on 
the  ones  that  yield  minimax  -  or  nearly  so  -  approximations. 

The  algorithms  are  presented  in  a  descriptive 
manner,  since  the  writer  had  to  rely  on  incomplete 
references  for  the  vast  majority  of  the  information 
contained  in  this  chapter.  It  was  not  easy  to  uncover 
the  assumptions  used  in  these  algorithms  or  to  evaluate 
their  efficiency  on  different  types  of  computers.  Thus 
consideration  of  these  points  has  been  omitted. 

The  problem  of  obtaining  a  minimax  rational 
approximation,  of  a  specified  pair  of  degrees,  to  a 
continuous  function  is  very  difficult,  because  it  is 
basically  a  non-linear  problem.  Some  algorithms  such 
as  the  extension  of  Remez's  Second  Algorithm  and  Maehly's 
Direct  Methods  attack  the  problem  directly.  But  other 
algorithms,  such  as  Maehly's  Indirect  Methods,  telescoping 
procedure  for  continued  fractions,  and  the  Linear 
Inequality  Algorithm  of  Loeb,  attack  the  problem  in  a 
manner  that  leads  to  either  a  correction  of  the  coeffic¬ 
ients  of  a  Pade  or  other  rational  approximant  to  give  a 
lower  order  approximant  that  is  nearly  a  minimax  one,  or 
transformation  of  the  non-linear  problem  into  a  linear 
problem.  Another  problem  -  that  of  assuring  an  approxima¬ 
tion  that  is  free  of  poles  in  the  interval  of  approximation 
-  is  overcome  in  many  algorithms  by  restricting  the 
denominator  to  values  that  are  greater  than  zero  for  any 
value  of  x  in  the  interval  of  interest. 


The  Fade  method  has  been  included,  because  Pade- 
approximants  are  often  used  as  the  initial  approximants 
of  Maehly's  telescoping  procedure  for  rational  functions 
and  Indirect  Methods. 
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Section  5-2  The  Pade  Approximants 

This  method  differs  from  the  one  given  in  Hall 
and  Knight  (see  page  369  of  [6])  in  that  it  allows  us  to 
obtain  rational  approximations  of  varying  pairs  of  degrees. 
On  the  other  hand,  the  method  given  in  Hall  and  Knight 
results  in  a  rational  approximation  in  a  restrictive 
continued  fraction  form  that  provides  convergents  in  which 
the  degree  of  the  numerator  is  alternately  one  greater 
than  and  equal  to  the  denominator. 

The  Pade  Table  is  a  general  method  by  which  any 
power  series  can  be  transformed  into  a  table  of  rational 
approximations  to  the  function  represented  by  the  power 
series . 


The  coefficients  p^.,  q^  of  a  rational  function 
pn(x)/qm(x),  where 


PnU)  = 


n 

2  p,  x 
k=o 


k 


and 

m  k 

q(x)  =  1+  2  q,x  , 

00 

are  deduced  from  the  c,  of  the  power  series  2  c,  x  , 

K  k=o  K 

for  the  function  to  be  approximated,  with  the  aid  of  the 


definition 
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(5.2.1) 


00  k 

q  (x)  2  c,  x  -  p  (x) 

mv  ’  k==0  k  *nv  ' 


xn+m+1  |  k 
k=o  k 


This  definition  yields  n+m+1  linear  equations  for  n+m+1 
unknowns  Pk,qj,(k  =  0,  1,  . n),  (j  =  1,  2,  .  m). 


since  q  =  1. 
o 


Pade  approximations  are  most  useful  if  n  =  m 
or  n  =  m+1. 

Taking  the  general  case  we  have  from  (5.2.1) 

0  ^  k  n  Yr  CO  vyi-j—  vo  — 1—  k-|-  "I 

.2  c.x^  2  q,  xK  -  2  p,  xK  =  2  7,  xm 

i=o  ^  k=o  ^  k=o  k  k=o  k 


Now, 


00  n  m  if 

2  c,x  2  q,x 
£=0  ^  k=o  ^ 


00  min_[f,m) 
=0 


2  * 


z,  q,  c„  ,  x 

l=o  k=o  ^  t-k 


so  that  we  get 


(5.2.2) 


00  min(ii,m) 
2  2 
$=o  k=o 


Qi  c,  ,  x 


J 


n  k 
2  Pi  x 


k  ,0-k  k=Q  ^  "  k=0  rk' 


1  Y  xm+n+k+l 


Now,  there  is  no  term  other  than  the  first  sum  on  the 
L.H.S.  of  (5.2.2)  that  has  powers  of  x  in  the  range 
n+1  <  &  <  n+m.  Hence,  replacing  i  by  n+s,  1  <  s  <  m. 
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we  have 


m 

(5-2-3)  klQ  Vn+s-k  -  °' 

which  we  can  write  as  a  system  of  equations  In  the 
following  way: 

(5.2.3a) 


cn 

cn-l 

*  ’  ‘  cn-m+l 

t 

i — 1 

a1 

L _ 

cn+l 

cn+l 

cn 

cn-m+2 

q2 

cn+2 

• 

cn+m-l 

cn+m-2 

...  c 

n  J 

.V 

-  cn+.m  - 

The  lowest  power  of  x  on  the  R.H.S.  of  (5.2.2)  is  n+m+1 . 
Hence,  we  have 


(5.2.4) 


k 

pk  =  ilo  +  ck-£ 


for  k  =  0,  1,  .  ..,  n.  The  highest  power  of  x  in  the 

second  sum  on  the  L.H.S.  of  (5.2.2)  is  n.  Hence,  letting 
l  =  n+m+k+l,  we  have 


n+m+k+l 

Tk  x 


v  «  v n+m+k+l 

i=o  ^i  cn+.m+k-i+l  x 


Hence, 


y 


k 


(5.2.5) 


m 

i=o  Cn+m+k-i+l  ’ 


.>!  lO  . 
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Thus  by  solving  the  system  of  equations  (5.2.3a)  for  the 
q^,  k  =  1,  2,  .  . . ,  m,  we  can  find  the  pk,  k  =  0,  1,  .  ..,  n, 

from  (5.2.4)  . 


The  y^.  provide  us  with  an  error  estimate. 
However,  a  crude  idea  of  the  error  may  be  found  from 
considering 


< 


2n+l 


For  this,  it  is  sufficient  to  compute  7  .  yQ  can  be 

evaluated  from  the  relation  y  =  A  /<5  ,  where 

o  n'  n 


A 


n 


C1  C2 


cn+l  cn+2 


'n+1 


'2n+l 


and  6n  is  the  principal  minor  of  An  obtained  by  omitting 

the  last  row  and  column  in  A  . 

n 

After  obtaining  the  p^  and  q^,  we  can,  if  we 
wish,  transform  our  rational  approximation  into  an 
equivalent  finite  continued  fraction. 


'  '  ;■  '  .  M 
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The  accuracy  of  this  type  of  rational  approxi¬ 
mation  decreases  rapidly  as  the  value  of  the  argument 
increases  because  the  expression  for  the  error,  |  En  |, 
usually  depends  much  more  on  the  smallness  of  the  reduced 
range  in  which  the  approximations  are  used  than  on  the 

factor  y  in  the  error  term.  Hence,  unless  we  have  a 
o 

strongly  convergent  power  series  it  is  not  advisable 
to  use  uncorrected  Pade  approximations.  We  will  examine 
methods  of  correcting  Pade  approximants  in  section  ^  A 
and  section  5 .6. 

Section  3.3.1  An  Extension  of  Remez^  Second  Algorithm 

This  algorithm  is  an  extension  of  Remez's 
algorithm  for  finding  polynomial  approximations  to  the 
determination  of  "best-fit"  rational  approximations. 

It  is  helpful  in  describing  the  algorithm  for 
minimax  rational  approximation  to  make  reference  to  the 
simpler  algorithm  for  polynomial  approximation  since 
many  of  the  details  are  exactly  the  same .  Defining 
r(x)  =  f(x)  -  pn(x),  a  <  x  <  b,  it  is  known  (see  Chapter 
II,  section  6)  that  corresponding  to  the  best  polynomial 
approximation  of  degree  n  to  f(x),  there  exists  a  set  of 
points  x^  and  a  minimum  deviation  p  such  that 
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(5-3.1) 


r(x±)  =  ( - 1 ) ±p  ,  i  =  0,  1,  ...,  n+1 


(5.3.2) 

|  r(x^)  |  >  |  r(x)  |  ,  a  .£  x  £  b  *  i  =  0,  1,  .  .  .,  n+1 


The  determination  of  the  set  of  points  x^,  the  .minimum 
deviation  p,  and  the  polynomial  p  (x)  of  best  approximation 
is  achieved  by  replacing  the  system  of  equations  (5-3.1) 
and  (5-3-2)  by  the  following  iterative  scheme,  in  which 

■fh  • 

the  superscript  j  applies  to  the  j  1  stage  of  the 
iteration . 

(5-3-3)  r'3'(xiJ)  =  (-1)1  pJ'  ,  i  =  0,  1,  .  ..,  n+1 

(5-3-^)  r^(x^~r''')  =  extremum  of  r^(x),  1  =  0,  1,...,  n+1 

and  for  at  least  one  i 

rJ(x^^°‘1)  =  max  |  r^(x)  | 

a<x<b 


The  solution  is  started  by  choosing  an  arbitrary  set  of 
points  (x^0)  and  using  the  resulting  linear  system  (5-3-3) 
to  determine  p°  and  pn(x) .  the  maxima  and  minima 
(including  the  greatest  extrema)  of  r°(x)  become  by 
(5-3-^)  the  set  of  points  [x^j  and  a  cycle  of  the 
computation  is  complete.  In  applications  the  process 


'  i'.j; sr;,:; 
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can  be  terminated  by  a  sufficiently  close  agreement  between 
successive  values  of  p  or  alternatively  when  there  is 
sufficiently  close  agreement  between  the  values  of  the 
|  r^(x^+1)  |  and  |  p1^  |. 


As  in  the  linear  case,  the  minimization  of 


P(Pn’qm> 


i  w  n  Pn(x) 

.max  fix)  -  7 — r 

a<x<b  ^m ( x ) 


can  be  achieved  by  the  minimization  of  an  auxiliary 
function 


max 

xeX 


I  f(x) 


PnU) 


for  an  appropriate  set  X  consisting  of  n+m+2  points. 

The  following  algorithm  systematizes  the  search  for  this 

set  X.  Given  a  trial  set  X,  select  p  °  and  q  °  to 

n  qn 

minimize  6.  If  6(p  °,q  °)  >  p(p  °,qm°)  then  (pn°,q  °) 
is  a  solution  to  the  problem  and  X  is  the  appropriate  set. 
Otherwise,  select  a  point  £  not  included  in  X  for  which 

1  f(i)  ■  1  =  A(Pn°'q'"°)  • 

Qm  '  ^  ' 


Next,  as  described  below,  replace  a  certain  element  x 
of  X  by  |  and  repeat  the  process  with  the  new  set  X. 

It  is  always  possible  to  do  this  in  such  a  way  that  the 
values  of  5  increase. 


. 

Jy.  .  :  >  o •  \-JC; ■  T,  \ 

Xiv’  '  ;/  :  -  '-±1  i.  v  :i,T 

T"  .0  .r;;. 

c  '■ '  'J o  si  sa  ,,  9a 
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We  can  modify  the  algorithm  In  the  following 

Vi 

way  to  give  faster  convergence.  At  the  r  stage  we 

r  r 

have  a  set  of  points  =  {x.,  ^..^^4.0}  that  is 


n 


n+m+2 
r 


used  to  generate  a  rational  function  Q^fx) •  Then 

T 

associate  with  each  point  x  .  that  extremum  in  its 

T 

neighborhood  nearest  to  x  .  at  which  the  error  has  the 

J 

P 

same  sign  as  at  Xj.  Denoting  this  set  of  extrema  by 


£t1'  Tn+m+23 


we  then  let  =  T^,  instead  of  replacing  one  xj  by 

that  member  of  Tr  at  which  the  error  has  its  greatest 
magnitude.  The  idea  behind  this  is  that,  since  the 
limiting  set  of  points  X  corresponds  to  the  extrema  of 
Q^^x)  (the  "best-fit"  rational  approximation),  convergence 
of  the  algorithm  should  be  more  rapid  using  the  above 
method  of  replacement. 


Actually,  we  do  not  need  a  rational  function 

Q^(x)  at  all  but  only  a  set  of  points  X1  and  a 

value  of  p  which  we  would  like  to  be  as  close  as 

* 


possible  to  X  and  Pnm^  respectively.  We  can  either 
use  for  X^  the  extrema  of  Tn+m+-^(x)  -  the  Chebyshev 
polynomial  of  degree  n+m+1,  or  the  extrema  of  an 
economized  approximation  to  f(x)  which  generally  also 
gives  a  good  initial  value  for  p . 


The  system  of  equations  (5.3.1)  is  non-linear 


in  the  rational  case.  If  we  write  them  in  the  form 


(5.3.5) 

[Ux^M-l)1?]  J1  qJ(x11)J  -Jo  PjU-lV  =  -f(x11)  +  (-l)± 


i 


0,  1 


•  •  •  ) 


n+m+1 


we  see  that  the  only  non-linear  factor  is  in  the  first 
term.  Now,  consider  the  system  (5*3.5)  split  into  two 
sets  of  equations,  one  of  which  we  denote  by  Sn+m+1 
consisting  of  n+m+1  equations  (which,  for  convenience 
is  taken  to  be  the  first  n+m+1  equations)  and  the  other 
being  the  remaining  last  equation.  The  system  Sn_j_m+^ 
is  linear  in  the  unknowns  A  =  {pQ, • . . * Pn, q^ j ^ Qm) 
and  may  be  thought  of  as  implicitly  defining  each  member 
of  A  as  a  function  of  p.  Therefore,  the  last  equation  may 
be  thought  of  as  a  function  of  p  alone  and  may  be 
written  as 


(5.3.6) 


F(p) 


11  1 

2  p  . (x  ,  , , 

j=o  n+m+1 


)J’  +  f(x. 


1 

n+m+1 


P 


0 


i  ■■  4  .  .j  ;  v 
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The  equation  F(p)  =  0  may  be  solved  by  the 
secant  method,  in  which  case  the  procedure  is  as  follows: 


(1)  Choose  two  values  of  p,  pQ  and  p-^  and  for  each  value 
solve  the  linear  system  Sn+m+1  for  pQ, ...,pn, 

^l*  *  ’  '  *  ’  Substituting  these  results  into  (5-3.6) 

we  get  F(  pQ)  anf  F(p  . 

(2)  Using  the  secant  method,  get  p2  as 


(5-3.7) 


=  p  + 


Pi  -  P 


o 


p(Pi) 


1  p(p0)  -  p(Pl) 

(3)  Solve  sn+m+i  with  p  =  p2  and  use  (5-3.6)  to  get 

f(p2). 

(4)  Repeat  steps  (2)  and  (3)  until  the  desired  degree 
of  convergence  is  achieved. 


In  practice  it  is  convenient  to  apply  (5.3-7) 
once  only.  This  can  be  done  because,  when  the  algorithm 
is  far  from  convergence,  high  accuracy  in  p  is  not  worth¬ 
while  and,  when  the  algorithm  is  near  convergence,  the 
values  of  p  change  very  little  from  one  stage  to  the  next 
which  enables  one  application  of  (5*3-7)  to  achieve  high 


accuracy. 
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The  key  to  the  remainder  of  the  computation  Is 

■jo 

the  search  for  Tr,  the  set  of  extrema  of  f (x)  -  Qnm(x) • 
Ralston  (see  [20])  has  used  a  search  prcedure  with  which 
he  has  searched  for  extrema  of  f(x)  -  Qnm(x)  by  starting 
at  x^r  and  searching  in  the.  direction  of  increasing 
magnitude  of  error.  When  the  peak  is  found,  he  then 
uses  quadratic  interpolation.  By  choosing  a  fairly 
coarse  step  size  at  the  first  stage  and  refining  this  as 
the  algorithm  converges,  the  total  computation  may  be 
kept  fairly  near  a  minimum. 

A  proof  of  the  convergence  of  this  algorithm 
has  been  given  by  Ralston  in  [20]  which  appears  to  place 
non-restrictive  conditions  on  f(x),  the  function  to  be 
approximated.  He  proves  the  convergence  of  the  algorithm 
by  assuming  the  existence  of  a  suitable  initial  approxi¬ 
mation  which  provides  the  first  extrema  of  r°(x)  with 
the  same  sign  as  the  first  extrema  of  r*(x). 

A  program  for  this  process  in  algorithmic 
notation  is  given  in  the  appendix. 

Section  5. 3. 2  The  Loeb  Weighted  Minimax  Algorithm 

We  wish  to  minimize  the  expression 

(5-3.8)  max  |  f (x)  -  Q(x)  |  , 

a<x<b 


where 
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n 


Q(x) 


Pn(x) 


qm(x) 


^  p  .X 

1=0  *n-l 

 ; 

i  +  ^2^  qm_i  x 


We  may  write  (5-3.8)  In  the  form 


max 

a<x<b 


qmM 


I  f(x)qm(x)  -  pn(x) 


since  qw(x)  must  not  vanish  in  [a,b]  and  hence  Q.m(x) 


xm 


can  always  be  made  greater  than  zero  throughout  [a,b] . 

We  can  now  use  the  following  iterative  procedure:  at 

the  k^h  step  select  pn^"(x)  and  qm^(x)  so  as  to  minimize 


max 


1  i  -o/ _ \~  k  /  —  \  -  k , 


a<x<b  q^-^(x) 
- Hm  v  ’ 


I  f(x)qmK(x)  -  pn  (x)  |  , 


n 


where  the  superscripts  k  indicate  the  stage  of  the  itera¬ 
tion.  In  this  subproblem,  l/[q^-1(x)]  is  regarded  as 
a  weight-function.  For  practical  purposes,  the  interval 
[a,b]  may  be  replaced  by  a  discrete  set  {x-j,  ,x^,  .  .  . , 


x 


3 


chosen  in  the  same  manner  as  the  initial  set 


n+m+2- 

of  ordinates  in  section  (5-3-1)  -  and  the  problem  is 
thereby  reduced  to  an  over  determined  system  of  linear 
equations  which  is  to  be  solved  in  the  minlmax  or  Chebyshev 
sense . 


The  advantage  of  this  method  is  the  linearity 


"■  '  o  •  .  "■EL  LiLt;.  Sfif 


'  :  nJ  ■  . 
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of  the  system  of  equations  which  arise  In  each  Iteration 
step . 


It  has  been  reported  by  E.  W.  Cheney  and 
T.  H.  Southard  (see  [4])  that  the  Space  Technology 
Laboratories  Inc.  has  used  this  algorithm  with  great 
success.  For  example,  it  is  reported  that  five  to  ten 
steps  of  the  iteration  easily  suffices  to  give  best 
rational  approximations  with  up  to  fifteen  variable 
coefficients . 

The  conditions  of  f(x)  which  are  sufficient 
to  quarantee  the  convergence  of  the  above  iterative 
process  are  not  known. 

Section  5-3.3  Loeb's  Linear  Inequality  Method 

This  method,  given  by  Loeb  in  [12],  permits 
us  to  use  linear  programming  techniques  to  obtain  a 
solution  to  the  problem  stated  in  section  5-3-2,  where 
the  expression  to  be  minimized  was 

p(pn,  qn,---,qm)  =  max  I  f(x)  -  Q(x) 

°  no  m  a<x<b 

Now,  if  the  number  p*  =  g.l.b.p  were  known,  the 
coefficients  Pj_,q_^  could  be  obtained  by  solving  the 


■ 
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the  system  of  Inequalities 

I  qm(x)f(x)  -  pn(x)  |  <  p*  |  qm(x)  | 

for  a<x<b .  Hence,  If  we  restrict  ourselves  to  functions 
for  which  qm(x)  >  a  >  0  in  [a,b],  we  may  write  the 
following  system  of  linear  inequalities 

qm(x)  >  a 

qm(x)f(x)  -  pn(x)  <  p*qm(x) 

Pn(x)  -  qm(x)f(x)  <  P*qm(x) 

for  a<x<b .  The  number  p*  is  not  usually  known  at  the 
beginning  and  has  to  be  determined  by  successive  trials. 

As  an  initial  value  of  p*,  say  p°,  we  could  take  the 
maximum  deviation  of  a  corrected  Pade  approximant  to  the 
function  in  the  interval  of  interest.  With  the  trial 
value,  p°,  the  systems  of  inequalities  may  be  tested 
for  consistency.  If  it  is  inconsistent  then  p°  <  p*; 
otherwise  p°  >  p* .  It  has  been  recommended  that  the 
consistency  test  be  made  by  minimizing  the  auxiliary 
function 

<5(p0*  •  •  .,qm)  =  max  max[a-qm(x),  |qm(x)f(x)-pn(x)  |-p°qm(x)  ] 
a<x<b 

The  system  of  linear  inequalities  above  is  equivalent  to 
the  single  inequality  6  <  0,  and  we  may  minimize  6  by 
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linear  programming. 

Section  5-3.^  A  Differential  Correction  Method 

This  method,,  which  was  given  by  Cheney  and  Loeb 
in  [3]*  is  effective  without  assumptions  on  f(x)  other 
than  continuity.  This  does  not  imply  that  all  of  the  other 
methods  discussed  require  more  restrictive  assumptions. 

At  the  start  of  the  process;  Q(x)  =  pn°(x)/qm°(x)  is 

arbitrary  except  for  the  condition  qm°(x)  >0  in  [a;b] . 
Now,  assuming  that  Q^(x)  has  been  determined;  put 

pk  =  max  |  f(x)-Qk(x)  |  . 
a<x<b 

We  then  determine  Q^+'L(x)  in  such  a  way  as  to  minimize 
the  expression 

6  =  max  [  |  f(x)q  (x)  -  p  (x)  |  -  p  qm  (x) ] 

a<x<b 

under  the  restriction  that  the  numbers 

I  pq  l  • .  -  .1  Pn  U  I,  • .  .,1  qm  I  be  <  1. 

It  has  been  shown  (see  [3])  that  i  g.l.b.p  =  p^ . 

The  calculations  to  determine  Q,k(x)  may  be  carried  out 
by  using  the  methods  of  "convex  programming"  since  5  is 
a  convex  function  of  the  coefficients;  and  the  constraint 


set  is  also  convex. 
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Section  5.4  Maehly*s  Telescoping  Procedure  for  Rational 


Approximations 

We  will  begin  by  examining  a  simple  economization 
process  of  a  polynomial,  and  then  examine  an  economization 
procedure  for  rational  approximants . 

Let  f(x)  be  an  analytic  function  which  we  wish 
to  approximate  on  an  interval  [x^,x2 ]  to  a  given  degree 
of  accuracy.  Let  us  assume  that  we  have  found  an  integer, 
n,  such  that  the  truncated  power  series. 


(5.4.1) 

n+1  ^ 

pn+l(x)  =  k^0  °kx  ’ 

is  a  sufficiently  good  approximation  on  our  interval.  We 
then  wish  to  find  a  shorter  polynomial  pn* (x)  so  that 
the  condition 


(5.4.2) 

max  |  p  *(x)  -  f ( x )  |  =  E  , 

X1<X<X 2  ^ 

where  E 

is  the  greatest  lower  bound  of  the  error  of 

the  approximation,  is  at  least  approximately  fulfilled. 
We  may  accomplish  this  by  the  "telescoping  procedure" 
which  can  best  be  described  if  we  write 


(5.4.3) 

X1  =  £U1  ;  X2  =  €U2 

where  u~ 

-  u-^  is  of  the  order  of  one,  while  e  is  a 
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parameter  which  characterizes  the  length  of  the  interval 
[x1,x2] . 

We  next  denote  by  ^n-(-]_(u)  the  polynomial  which 
is  uniquely  defined  by 

(5.4.4)  S  ,  ,  (u)  =  2  S,  uk  +  un+1 

'  n+lv  '  k=o  k 

and 

(5.4.5)  max  |  S  +1(u)  |  =  Eg  , 

ui<u<u2 

where  Eg  is  the  greatest  of  the  deviations  of  ^n+]_(u) 
from  zero  in  [u^u^]  .  Sn+1(u)  essentially  a  Chebyshev 

polynomial . 


The  telescoping  procedure  consists  of  taking 
Pn*U)  =  Pn+1(x)  "  cn_j_i 1 SR+ 1 ( u )  where  u  =  x/e  .  If 
the  power  series  converges  for  x^<x<x,~,  3  ^hen 

(5.4.6) 

pn*(x)-f(x)  =  Cpn*(x)-Pn+1(x) ]  -  [f (x)-Pn+1(x) ] 


- c  ,  ,  en+^S  , (u)  - 
n+1  n+lv  ' 


_n+2 


00  k 

2  e  c  .  0  . ,  u 
k=o  n+2+k 


n+2+k 


from  which  it  follows  that 


(5.4.7) 


lim 

e-»o 


pn*(eu)-f (eu) 
en+l 


"cn+lSn+l 


(u) 


■x-' v r';  /  :!*  .  r-i:\ 


;  -  ■ ; :  ■'  '  l '  |  '•  '  a*  '■  - xsn  ,  . 1 


,  • 


■  a 3. 
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The  theory  of  Chebyshev  polynomials  tells  us 
that  Sn+i(u)  assumes  the  maximum  of  Its  absolute  value, 
with  alternating  signs,  exactly  n+2  times  within  the 
interval  [u,,Ug]  .  Now  (5-4.7)  shows  that  this  is  also 
at  least  approximately  true  for  pn*(eu)  =  f(eu)  if 
e  is  sufficiently  small;  however  this  property  characterizes 
the  best-fit  polynomial  of  degree  n  of  our  function  f(x) 
uniquely,  according  to  a  fundamental  theorem  of  Chebyshev 
(see  chapter  II,  section  6).  Hence  it  is  fair  to  say 
that  (5.4.2)  is  at  least  approximately  fulfilled  if  cn+-^  ^  0 
and  if  e  is  sufficiently  small. 

The  telescoping  procedure  can  be  applied  almost 
as  easily  to  any  type  and  form  of  continued  fraction. 

Suppose  that  the  function  f(x)  which  we  wish  to  approxi¬ 
mate  can  be  represented  by  a  convergent  continued  fraction 
of  the  form 


(5.4.8) 


f(x) 


qp  |  +  aix  |  +  a2x  | 

I  ^O  I  b^  i  bg 


+  . .  .  + 


ak*  1 


+  .  .  . 


and  let  us  assume  that  the  (n+l)th  approximant,  which  can 
be  written  as 


R 


n+1 


(x)  = 


a 


+ 


OtjX  |  OgX 


+ 


o 


b- 


b, 


+ 


,  anx  I  ,  an+lx 
+  -  + 


n 


b 


n+1 


(5.4.9) 
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represents  f(x)  well  within  the  given  error  bounds  on  the 
interval  [x^x^]  .  Our  wish  is  to  find  a  shorter  approxi- 
mant,  Rn*(x),  so  that  the  condition 

(5.4.10)  max  |  R  *(x)-f(x)  |  =  E  , 

x-^<x<X2 


where  ER  is  the  greatest  lower  bound  of  the  error  of  the 
approximation,  is  at  least  approximately  fulfilled. 


Let  us  first  look  at  the  uncorrected  n"^h 
approximant  Rn(x) .  From  the  theory  of  continued  fractions 
it  is  known  that 


(5.4.11) 


lim 

e-*o 


Rn(eu)  -  f(eu) 
€n+l 


cn+lx 


n+1 


where  cn+]_  is  given  below  by  (5.4.16).  Now  (5.4.11)  is 
fully  analogous  to  (5-4.7).  Hence  we  can  hope  to  find 


corrections  to  the  coefficients  of  Rn(x)  such  that 


* 


(5.4.12) 


lim 

e->o 


Rn  ( eu)  -  f(eu) 


.n+1 


cn+lSn+l 


(u)  , 


where  u,  e,  and  Sn+^(u)  are  as  defined  above. 

Half  of  the  constants  in  (5.4.8)  and  (5-4.9) 
are  redundant;  two  continued  fractions  which  have  the 
same  quotients  r^.. 


v-  '!-  |  j  y  ■■■  ,}  y :>y.  <7 

•  . 
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(5.4.13) 

rt  -  “f1  ,  k>0  , 

k  bkbk+l 

and  the  same  cxq,  are  exactly  equivalent.  Hence  to 
* .  . 

obtain  Rn  (x)  we  may  alter  either  the  a  of  the  b^  of 

R  (x) .  In  the  first  case  we  take 
n\  ) 


(5.4.14a) 

a  *  I  a  *x  I  a  *x  I 

R  a  (x)  =  0'+1  +  ...  +  n 

n  1  bo  1  bl  1  bn 

(5.4. l4b) 

V  “  ak  +  V  ’ 

(5 .4. l4c) 

„  1  q  /  _  \  n+ 1  -  k  n 

k  “  kSk  ^  ^  n,  rm  • 

m=k 

The  corresponding  formulae  for  correcting  the  b^  are: 


(5.4.15a) 

D  b*,  x  ao  i  .  alx  I  ,  ,  anx  1 

n  b  *  b*  b  * 

o  l  n 

(5.4.15b) 

b  *  =  b,  +  b,  ' 

k  k  k 

(5.4.15c) 

V  =  bksk  (-)n+1-k  n  rm 

m=k 

For  finite  values  of  e  the  two  forms  yield  slightly 
different  numerical  values,  but  both  of  them  satisfy 
(5- 4.12)  in  the  limit  e-*o .  It  may  also  be  shown  that 


.  ;  , . . 


' 


.  :  ...  ,  (/. 
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(5.4.16) 


n+1 


a  n 

(_l)n  o  n  r 

v  '  bQ  m=l  m 


An  approximant  of  the  form  (5.4.8)  can  also  be  expressed 
as  the  quotient  of  two  polynomials. 


(5.4.17) 


Vx) 


pnW 


QLn(x) 


which  can  be  computed  from  the  recurrence  formulae 


(5.4.18) 


p  =  bp  .,+ap 


m 


rrrm-1  nrm-2 


q  =  bq  n+aq  _ 
Hm  mHm-l  mHm-2 


starting  with 


for  m>2 


(5.4.19) 


P  =  cx  ,  p,  =  a  b, 
*o  o’  ^1  o  1 


qO  =  V  ql  =  Vl  +  alx 


The  rational  approximation  satisfying  (5.4.12)  can  be 
found  from 

y* 


n 


(5.4.20a)  Q^hx)  = 


P 


n 


q. 


T 

n 


Pn+Vx  kg2  YkPk-2 

%+rj^x  +  x^S  'Vicqk-2 


with 


J  n  1  n+l  a 

-  -sk(-€)  \  £ 

K  K  m=k 


(5.4.20b) 


. 

■  '  1  .  ;k  '  . 
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The  last  form,  (5.4.20),  Is  probably  the  most  convenient 
for  practical  purposes,  especially  if  the  division  time 
-  of  the  computer  -  is  more  than  twice  the  multiplication 
time  . 


Ralston,  in  [19],  has  given  another  method  for 
the  economization  of  rational  approximations  that,  unlike 
Maehly's  procedure,  can  be  applied  to  all  Pade  approxima¬ 
tions  .  The  application  -  in  some  cases  -  of  the  process 
is  simpler  than  that  of  Maehly’s  procedure,  because 
Ralston  uses  the  coefficients  of  the  Chebyshev  polynomial 
(of  the  first  kind)  to  obtain  the  corrections  to  the 
coefficients  of  the  approximation  that  give  a  minimax  - 
or  nearly  so  -  rational  approximation. 

Section  5-5  Maehly's  "Direct  Methods"  for  Fitting 

Rational  Approximations 

The  term  "direct  methods"  indicates  that  the 
coefficients  of  the  Chebyshev  -  approximant  are  computed 
directly.  Maehly  uses  a  weight  function,  g(x),  and  thus 
his  condition  for  the  "best-fit"  Chebyshev  -  approximant 
to  a  function,  f(x),  is  that 

|  Q^x)  -  f(x)  | 
g(x) 


(5.5.1) 


max 

a<x<b 


g(x)  >  0, 


-  ; 
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OI;::r  Hi  X  C.X  l  l  XXl  -3:1.;  ,  •  , 

.'1  ;  '  •  ■  XO'l' "\0  IX.XS.:  0  I  V: 

i 

■  '  .  .  ..  '  .''2.2  ..  ■  ■  .  J.  y  :  r  1  :  ■  -  • 


■  '  .  .  ■'  V." . ■■  i  :  :  '■  '  ■'  j.:  ,:.i  „  ^ .{  ;t ...  :  .3.  j . 


’  .  .  noi:  1?.  ' 


131. 


be  minimized,  where  (x)  is  the  set  of  rational 

functions  of  the  form 


(5.5.2) 


Pn(x)  _  Pn+Pn-lx+- • •+poxn 
%,(x)  ”  •  -+<:loxra 


Each  rational  function  (x)  gives  rise  to 

an  error  curve,  or  error  function,  defined  by 


(5.5.3) 


r(x) 


-  f(*> 
gU) 


and  we  write  the  error  curve  of  the  "best-fit"  approximant 
as 

(x)  -  f(x) 

g(x) 


(5.5.4) 


r* (x)  = 


* 


m 


An  error  curve  is  said  to  have  standard  form 
if  it  satisfies  the  following  additional  conditions: 


(1)  there  are  exactly  n+m+2  critical  points, 

a=x*<x*<  ...  <  x*  =  b  • 

o  i  n+m+1 

(2)  |r(x)  |  =  p  at  the  critical  points  with  alternating 

signs . 

(3)  r*(x)  is  continuous,  and  vanishes  only  for 


=  x. 


X 


i  =  1,  2,  .  .  ., 


n+m . 


. 


‘1 


! 
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A  sufficient,  but  not  necessary,  condition  that  the 
error  curve  r(x)  be  optimal.  Is  that  it  has  standard 
form.  All  the  methods  discussed  in  this  section  are 
based  on  the  assumption  that  r*(x)  has  standard  form. 


The  unique  solution  r(x)  =  r*(x),  x^  =  x^' 
of  the  following  set  of  equations 


(5.5.5)  r(xjL)  =  (-l)p,  i  =  0,  ...  ,  n+m+1 


d 


^  r(xi)  =  0,  i  =  1,  ...  ,  n+m 


or 


(5-5.6)  Qnm(x1)  -  f(xj_)  -  (-l)1pg(x1)  =  0, 

1=0,  ...  ,  n+m+1 

(5.5.7) 


g(xl>  ((Wxi>-f(x:L>>  -  (+m(xh-f(xi))  ai  s(xP  =  °. 


i  =  1,  ...  ,  n+m, 

is  supplied  by  every  standard  error  curve  and  its  critical 
points.  There  are  2(n+m+l)  equations  and  the  same  number 
of  unknowns.  Hence  the  problem  consists  in  solving  the 
above  system  of  nonlinear  equations . 
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I .  The  First  Direct  Method 

This  is  a  two-stage  iteration  method  which  can 
be  described  as  follows: 

(1)  Make  an  initial  guess  at  the  critical  points  and 
solve  the  equation  (5.5*6)  to  obtain  a  rational 
function  Qnm(x)  an^  a  value  of  p. 

(2)  A  new  guess  at  the  critical  points  is  then  given  by 
the  extrema  of  r(x),  and  this  guess  is  used  in  the 
subsequent  stage  (l)  step.  It  also  solves  the 
equations  (5.5*7). 

For  stage  (2)  a  simple  searching  procedure  is 
recommended  which  is  also  used  in  the  Second  Direct 
Method  and  will  be  described  later. 

The  determination  of  the  polynomial  equation 
for  p  poses  problems,  expecially  for  large  m;  also, 
the  evaluation  of  the  equation  is,  at  times,  difficult. 
It  is  also  at  times  difficult  to  determine  the  value  of 
p  that  provides  a  pole-free  approximant .  Hence  use 
of  the  First  Direct  Method  is  not  recommended. 

II .  The  Transformation  of  Rational  into  Polynomial 
Approximation 


In  the  polynomial  case,  stage  (l)  of  the  First 
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Direct  Method  is  greatly  simplified  because  the  equation 
for  p  is  linear.  Hence,  it  is  useful  to  note  that 
fitting  a  rational  approximant  (x)  can  1°e  reduced 

to  fitting  a  sequence  of  polynomials  of  the  same  degree. 


Choose  two  reference  polynomials  p  (x)  and 
q  (x)  such  that  the  quotient  pn(x)/qm(x)  is  irreducible  and 
at  least  one  of  them  actually  attains  the  degree  n  or  m, 
respectively.  The  optimal  error  curve  may  then  be 
written  as 


(5.5.8) 


*/  \  1 
r  (x)  = 


Pn(x)  pn(x) 


n 


g(x)  q*(x)  q  (x) 


)  ~  (f(x)  - 


Pn(x) 


q  (x) 


■)] 


(p*(x)qm(x)“q*(x)pn(x) )  -  U)qm(x)-pn(x) ) 


g(x) q  (x ) q  (x) 
m v  '  ^m v  ' 


where  Q^m(x)  =  Pn(x)/o.m(x)  denotes  the  Chebyshev 
approximant  of  f (x)  with  the  weight  function  g(x) . 
Clearly  we  may  say  that 


(5.5.9) 


»n+m(*>  -  F(x) 

H(x) 


r*(x) 


is  also  the  error  curve  of  the  polynomial 

f  i  (x)  =  p*(x)q  (x)  -  q*(x)p  (x) 

rn+mv  '  ^nv  '  nrr  '  vrr  /J"nv  ' 


(5.5.10) 


: 
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with  respect  to  the  function 

(5.5.11)  F(x)  =  q*(x) (f (x)qm(x)  -  Pn(x)) 

and  the  weight  function 

(5.5.12)  H(x)  =  g(x)q*(x)qm(x) . 

^n+m^x^  "kest-fit"  polynomial  since  r*(x)  has 

standard  form.  Hence  we  may  attack  our  problem  by 
Remez  *s  Second  Algorithm.,  for  example. 

*  .  . 

Since  q  (x)  is  not  known  at  the  beginning  of 
the  iteration,  we  must  start  with  a  guess  for  q*(x) • 

Then  one  step  is  taken  in  the  direction  of  the  "best-fit" 

*  /  X 

polynomial  fn+m\x) •  From  the  linear  system 

<5-5-13)  K+mW  =  PnU)qm(X)  -  qmMPnU) 

we  determine  the  new  polynomials  Pn(x)  an(l  qm(x), 
which  are  approximations  to  p*(x)  and  q*(x),  respec¬ 
tively.  The  new  q  (x)  enters  the  subsequent  iteration 
step.  The  linear  system  (5. 4.13)  has  n+m+1  equations  and 
n+m+2  unknowns,  but  the  quotient  pn(x)/qm(x)  has  in 
fact  only  n+m+1  free  coefficients.  Hence  we  may  fix  one 
of  the  coefficients  of  Pn(x)  or  qn(x)  and  then  the 


«  '  v'?c  ;  ' 
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system  can  be  solved.  It  Is  Important  that  the  reference 
polynomial  q  (x)  be  positive  In  the  Interval  [a,b], 
and  p  (x)/q  (x)  should  not  be  a  good  approximation  to 
f(x);  for.  If  It  Is,  the  formation  of  the  difference 
f (x)qm(x)-pn(x)  may  result  in  a  substantial  loss  of 
accuracy. 

Ill .  The  Second  Direct  Method 

The  error  curve,  r*(x),  has  exactly  n+m+1 
zeros,  z1*,  1=0,  ...,  n+m  in  the  interval  [a,b]  if 

it  is  of  standard  form.  Hence  we  may  write 

n+m  * 

(5.5.14)  r*(x)  =  G(x)  n  (x  -  z  ) 

k=o  K 

If  we  characterize  r(x)  by  its  zeros  rather  than  by 
its  extrema,  we  avoid  the  instability  of  the  First 
Direct  Method.  This  leads  to  the  Second  Direct  Method 
which  can  be  described  as  follows: 

(l)  Letting  a<z  <  ...  <  z  <bbea  guess  at  the 

o  n+m 

*  .  . 

zeros  of  the  optimal  error  curve  r  (x),  we  obtain, 
through  rational  interpolation,  an  (x)  which 

satisfies  the  conditions 

=  f ( zk)  j  k  =  0,  . . . ,  n+m. 


(5.5.15) 


zk^ 


:  ,::;og  •  .rIT  :;:a 
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(2)  The  zeros  of  r(x)  are  corrected  by  the  calculated 
extrema.  The  corrected  zeros  k  =  0,  .  ..,  n+m 

are  then  used  in  the  subsequent  stage  (l)  step. 


The  main  problem  of  this  method  Is  the 

correction  of  the  zeros  z^  using  the  extrema  x^. 

Maehly  used  the  following  method  for  correcting  the 

zeros  z^  which  he  obtained  through  use  of  a  short 

variational  argument  (see  page  265  of  [ 16 ] ) .  The 

correction.  6z,  .  to  the  zeros  z,  can  be  found  from 
’  M.  k 

the  following  system  of  equations 


n 

2 

k=o 


xo"zk 


)  6z 


k 


Sin 


rUi) 

r(xo) 


i  =  1,  .  .  .  ^n-fm+1^ 


where  we  may  replace 


in  |  - |  by  2  J — — i - 1 — -11— 1 

r(xQ)  I  r(xi)  I  +  I  r(xQ)  I 

Hence  our  system  of  equation  becomes 

n  (xQ-x1)6zk  =  2  I  WXj)  I  -  I  r(x0)  I 

k=°  (x1-zk)(xo-zk)  I  rUp  |  +  |  r(xQ)  | 


ly  ...  y  n+rrri-1 


which  is  well-conditioned  since  it  has  its  largest 
elements  close  to  the  diagonal. 
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In  the  determination  of  the  extrema  of  the  error 
curve  r(x)  we  cannot  use  Newton's  method  because  the 
error  curve  inevitably  carries  "noise".,  which  precludes 
the  numerical  computation  of  the  derivatives  r'(x)  and 
r"(x).  Hence  a  searching  procedure  was  adopted  which 
searched  in  equal  distances  for  the  largest  or  smallest 
value.  As  soon  as  a  value  r(x)  was  found  which 
surpassed  its  neighbors  r(x-h),  r(x+h)  the  searching 
was  stopped  and  the  three  points  were  interpolated  by 
a  parabola  with  extremum  x  where 


'■v 

X 


x 


^r (x+h) 

r (x+h) 


r (x-h) 
2r(x)  + 


_ )  h 

r(5c-h)  ^ 


The  selection  of  the  searching  distance  h  must  be 
given  considerable  care  since  the  bulk  of  the  computation 
in  the  direct  methods  consists  of  evaluating  the  function 
f(x)  . 


Section  5.6  Maehly's  "Indirect"  and  "Combined"  Methods 

for  Fitting  Rational  Approximations. 

The  term  "indirect  methods"  indicates  that  the 
coefficients  of  a  given  approximant  are  corrected  by  the 
addition  of  suitable  quantities  to  give  near  "best-fit" 
approximants .  The  indirect  methods  described  in  this 
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section  require  the  representation  of  a  function  f(x), 
by  a  finite  continued  fraction.,  to  the  full  accuracy 
required : 


(5.6.1) 


f(x) 


+ 


«  lx  I 


+ 


+ 


The  function  is  to  be  approximated  by  a  rational  function 
Rv  (x),  with  v  <  N,  in  the  Chebyshev  sense.  This,  like 

the  telescoping  procedure  for  section  5  A,  can  be 
regarded  as  a  problem  of  economizing  continued  fractions. 


The  original  Indirect  Method  dealt  with  the 
case  N  =  v  +  1  but  we  shall  discuss  the  Combined 
Methods,  which  deal  with  the  general  case  N  >  v.  Let 
us  first  describe  the  Combined  Methods  for  polynomial 
approximation  which  are  simple  and  may  serve  as  an 
introduction  to  the  more  complicated  rational  function 
case . 


I .  The  Combined  Methods  for  Polynomial  Approximation 

Assuming  that  f(x)  is  given  as  a  high  degree 
polynomial. 


N 

2 

k=o 


k 


ckx 


(5.6.2) 


f(x) 


within  the  interval  [o,b] ,  we  look  for  the  polynomial 
Pv  (x)  which  approximates  the  function  f(x)  best  in 
the  Chebyshev  sense . 


The  nth  Pade-polynomial  p,,(x)  =  2  c,  xk 

v  k=o  k 

is  a  good  approximant  for  f(x)  if  b  is  sufficiently 
small.  py  (x)  is  characterized  by  the  error  curve 


(5.6.3) 


r*(x) 


Pv*(x)  -  f ( x ) 

g(x) 


having  standard  form.  If  we  write 


Ap v ( x )  =  Pv*(x)  -  pv(x) 

we  may  write  (5.6.3)  in  the  form 

*  Ap v ( x )  -  (f(x)-pv(x)) 

r  (x)  =  - 

g(x) 

which  shows  that  Apv(x)  is  the  Chebyshev  -  approximant 
of  f (x)  -  pv(x)  with  respect  to  the  weight  function  g(x) . 
A  crucial  point  is  to  find  a  method  to  evaluate  the 
difference  F(x)  =  f(x)  -  pv(x)  without  a  substantial 
loss  of  accuracy. 


'  5U uYi -V1'  '' 


' 


l4l . 


Obviously,  in  the  polynomial  case,  we  have 

N  k 

P(x)  =  S  c,xK  . 

k= v+1  K 

Hence  any  direct  method  can  be  employed  to 
determine  a  "best-fit"  polynomial  Apv(x)  to  F(x), 
which  is  then  added  to  the  Pade-approximant  pv(x)  . 

II.  A  Formula  for  f(x)  -  Qv(x) 

If  f(x)  can  be  represented  in  the  form 


f(x) 


+  . . .  + 


the  combined  methods  consist  in  splitting  off  the  n^ 
Pade-approximant 


Qv(x)  = 


a 


o 


b 


ax  | 

+  — —  + 
lbl 


+ 


avx 
I  bv 


The  crucial  point  is,  as  above,  to  find  a  formula  for 
Qv(x)  -  f(x)  which  can  be  evaluated  without  a  substantial 
loss  of  accuracy. 


Referring  to  the  formulae  derived  in  section 
4.2  we  may  express  our  continued  fraction  as  convergents, 
p^/q^3  which  are  obtained  from  the  recursion  relations 
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pi  aipi-2  +  bipl-l 
qi  =  ai^i_2  +  biqi-l 


starting  with 


P_2  =  1*  P_x  =  0,  q_0  =  0,  q_1  =  1  . 


-2 


L-1 


Maehly,  in  the  following  manner,  derived  an  expression 
for  the  difference  Py/PjNj  ~  Pv/qy  *  N  >  v  .  Putting 


/  \  _  av+lx  *  ,  av+2x 

v+l,N(x)  |by+1  +  |  b. 


+  .  .  .  + 


av+NX 


v+2 


b 


v+N 


we  may  write 


(5.6.4) 


'N 


a 


o 


q. 


N 


b 


+ 


axx 


b 


+ 


+ 


1 


avx 

I  b.. 


+ 


v+l,N 


x) 


1 


An  expression  for  the  difference  of  two  successive 
convergents  can  be  written  as  (see  section  4.5) 


pi+i  _  £i 
qi+l  qi 


p1+lqi  -  p.q.+1 
qi+lqi 


qi+iqi  k=° 


1+1 

n  (-akx)  = 


bi+l 

‘i'-i-i  a1+1 


q, (q,  n+ 


q1)  k=o 


1+1 

n  ( _akx) 
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Applying  this  formula  to  (5-6.4)  we  obtain 


Qi 


vVHv-l 


+ 


v+l,N 


(x) 


V 

n 

k=o 


(-vo  • 


Hence 9  if  we  let 


sn(x)  =  4„.i  and  qm(x)  =  q„  , 


‘v-l 


m 


we  obtain 


(5.6.5) 


f(x)  -  Qv(x) 


(-x) 


v+1 


qm(x)(Sn(x)  +  _bU) 


V 

n  a 
k=o 


k  ' 


n 


v+l,N 


(x) 


which  is  the  desired  expression. 


Ill .  The  Combined  Method  for  Rational  Approximations 

The  combined  method  described  here  arises 
from  the  First  Direct  Method.  Essentially  the  method 
amounts  to  choosing  polynomials  p  (x)  and  qm(x), 
whose  quotient  is  the  Pade-approximant  Q  (x),  as 
reference  polynomials,,  and  proceding  as  in  section  5. 5.  II. 
We  recall  that  the  error  curve, 

Q*(x)  ~  f(x) 

g(x) 


r*  (x) 


■  -  ■;  ^  \  t  : 
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was  conceived.  In  that  case,  as  the  error  curve  of  the 
polynomial 

(5.6.6)  V'yW  =  pn(x)%i(x)  "  qm(x)pn(x) 
approximating 

F(x)  =  q*(x)  (f (x)qm(x)-pn(x) ) 

with 

H(x)  =  g(x)qm(x)qm(x) 

as  the  weight  function. 

Contrary  to  the  direct  approach  of  section 
5. 5- II.,  we  insist  that  pn(x)/qm(x)  be  a  good  approxima¬ 
tion  to  f (x) .  This  is  possible  because  formula  (5-6.5) 
enables  us  to  avoid  seperate  evaluation  of  f(x)  and 
Q  (x) .  If  we  substitute  the  correction  polynomials 

Ap  (x)  =  p  (x)  -  p  (x),  Aq  (x)  =  q  (x)  -  q  (x), 

into  (5.6.6)  we  get 

(5.6.7)  =  Apn(x)qm(x)  -  Aqm(x)pn(x) 

from  which  Apn(x)  and  Aq^(x)  can  be  computed  directly. 
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Since  the  equation  (5.6.7)  allows  one  degree 

"X-  *x* 

of  freedom,  we  normalize  p  (x)  and  q  (x)  by  reQuirinS 

■jf .  .  ^ 

that  q  (x)  have  the  same  constant  term  as  the  Pade- 
polynomial  q  (x),  which  leads  to  the  additional  relation 

(5.6.8)  Aqm(0)  =  0. 

The  method  is  again  a  two  stage  algorithm 
and  can  be  described  as  follows: 

(l)  A  guess  at  the  critical  points  is  made  for  which  we 
determine 

Qv(x±)  =  Pn^fAlm^P 

such  that 

Qv(x±)  =  (f (xi)-Qv(xi))  +  (-1)1  pg(x.)  , 

i  =  0,  .  .  . , n+m+1 . 

The  iteration  runs  as  follows:  make  an  initial  guess 
at  qm(x)  an-h  compute  f v(x)  and  p  from 

(5.6.9)  f  v(x^)  =  P(x±)  +  (-l)1pH(xi)  . 

We  then  determine  Ap  (x)  and  Aq  (x)  such  that 

tpv(x)  =  ApR(x)qm(x)  -  Aqm(x)pn(x)  . 


(5-6.10) 


■  l  '  '  l  ■  :  ; 
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Finally,  we  replace  q  (x)  by  q  (x)  +  Agm(x)  and 
return  to  the  beginning  of  the  iteration. 

(2)  The  extrema  of 

.  .  ~  F(x) 


are  used  as  a  new  guess  at  the  critical  points  and  stage 
(l)  is  re-entered. 

As  an  initial  guess  at  the  critical  points 
we  use  the  critical  points, 

x ±  =  (1  -  cos  |  ,  i  =  0,  ...,  n+m+1, 

of  the  transformed  Chebyshev  polynomial,  if  the  interval 
[0,b]  is  not  too  large.  The  initial  guess  at  q  (x) 
is  taken  to  be  q  (x),  the  denominator  of  the  Pade- 
approximant  Qy(x) . 

Maehly  suggests  the  following  iteration 
pattern  since  experience  has  shown  that  the  initial 
guess  at  q*(x)  is  m^ch  poorer  than  the  initial  guess 
at  the  critical  points.  Starting  with  two  or  three 
steps  of  the  stage  (l)  iteration,  keeping  the  initial 
x^,  we  obtain  a  Qm(x)  close  enough  to  Q^(x)  so  the 
next  stage  (l)  step  may  be  combined  with  one  step  of  the 
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stage  (2)  iteration.  Thus  the  combined  iteration  step 
consists  of  determining  fv(x)  by  (5*6.9)*  using  the 
resulting  error  curve  for  correction  of  the  x^  in 
stage  (2),  and  finally  correcting  q  (x)  . 

The  method  used  for  the  determination  of 
^  (x)  by  (5.6.9)  was  the  following  modification  of 
Newton's  interpolation  method.  Consider  two  polynomials 
cv+1(x)  and  Dv+1(x)  which  satisfy 

(5.6.11)  cv+l(xi)  “  F(xf  h 

(  1=0,..., n+m+1 

Dv+l(xf  =  (-D1H(x1)) 

We  then  have 

lMx)  =  cv+l(x)  +  pDv+l(x) 

where  p  is  uniquely  determined  because  the  highest 
powers  of  Cv+1(x)  and  Dv+1(x)  must  cancel. 

The  two  linear  systems  (5. 6. 11)  can  be  solved 
simultaneously  for  the  coefficients  of  Cv+^(x)  and 
Dv+]_(x)  since  the  systems  differ  only  with  respect  to 
their  right-hand  sides.  The  polynomials  are  written  in 
Newton  form. 
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v 


<WX)  =  co  +  cl(x'xo)  +---+  cv+l  150(x“xi) 


V 


3v+1(x)  =  an  +  di(x-xn)  +•••+  du+1  ,5„(x-xi)  * 


o  1 v  o 


v+1  i=ov  i' 


in  order  to  improve  the  condition  of  the  linear  systems. 
From  this  we  obtain  a  triangular  linear  system  with  two 
right-hand  sides  for  the  coefficients  c.  and  di . 

We  also  have  p  =  “cv+]/dv+]_  •  0ur  resulting  polynomial 

^v(x)  again  has  Newton  form, 

v- 1 

■fvW  =  eo  +  efx-x0)  +---+  ev  i£o(x-xf  ’ 


which  guarantees  stable  evaluation  of  the  error  curve. 
However,  ^v(x)  must  be  converted  into  the  customary 
polynomial  form  for  the  computation  of  the  corrections 
Apn(x),  Aqm(x). 

Let  us  now  consider  a  linear  system  for  the 
correction  polynomials.  The  coefficients  of  the  reference 
polynomials  Pn(x)  an<^  c4n(x)-’  (^e'i:ermine  the  coefficients 
of  the  linear  system  (5.6.10).  This  can  pose  serious 
problems,  in  most  cases,  because  the  coefficients  of  the 
reference  polynomials  display  different  orders  of  magnitude. 
Maehly  solved  the  problem  by  triangularizing  the  system 
(5*6.10)  in  such  a  way  that  it  can  be  solved  without  a 
substantial  loss  of  accuracy.  For  a  detailed  description 
of  the  method  used  see  pages  272-273  of  [ 16] . 


qA 


:  :  f-r  < 

"■  lo  'dr-.  '  .i.,,.  ■  '  ,;t  ■  ,  ^  /p 


149. 


CHAPTER  VI 

CONCLUSION 

In  this  thesis  an  attempt  has  been  made  to 
employ  basic  Chebyshev  concepts  In  the  general  case  of 
approximation  by  rational  functions. 

It  was  seen  in  chapter  III  that  the  process 
of  fitting  a  rational  approximation,  say  pn(x)/qm(x), 
to  a  prescribed  set  of  n+m+2  points  gives  rise  to  m+1 
values  of  minp .  At  most  one  of  the  values  of  mlnp 
will  give  us  an  approximation  that  is  free  of  poles  in 
the  interval  of  interest.  A  criterion  that  tells  us, 
when  we  have  a  set  of  four  ordinates,  whether  or  not 
we  can  obtain  a  simple  rational  function  of  the  form 

a0x  +  a]_ 

b  x  +  b.  ' 
o  1 

corresponding  to  either  value  of  minp,  that  has  no 
pole  in  the  interval  of  approximation,  was  given  -  in 
the  form  of  a  theorem  -  in  chapter  III.  Prom  the  writer's 
experience  with  higher  order  rational  approximations, 
it  appears  that  the  criterion  can  be  generalized. 
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In  chapter  IV  It  was  seen  that  It  is  not  always 
possible  to  obtain  an  acceptable  rational  function  that 
passes  through  a  prescribed  set  of  n+m+1  points  and  is 
of  a  specified  form;  that  is,  the  function  may  have  a 
pole  in  the  interval  of  interest,  or  the  system  of  equations 
for  obtaining  the  coefficients  of  the  rational  function 
may  be  inconsistent.  The  fact  that  the  system  of  equations 
may  be  inconsistent  can  raise  difficulties  in  some  of  the 
algorithms  of  chapter  V. 

Of  the  various  algorithms  for  obtaining  inter- 
polatory  rational  functions,  Thacher  and  Tukey's  method 
seems  preferable  to  the  others.  The  two  main  reasons 
for  this  preference  are  the  generality  of  the  algorithm 
and  the  fact  that  the  algorithm  does  not  require  a  large 
amount  of  storage  for  intermediate  results  when  it  is 
programmed  for  a  computer. 

The  algorithms  for  obtaining  minimax  rational 
approximations,  which  are  described  in  chapter  V,  can 
be  considered  as  belonging  to  two  classes;  that  is, 
methods  that  take  the  direct  approach,  and  those  that 
take  the  indirect  approach.  It  appears  that  use  of  the 


methods  of  the  first  class  will 


in  most  cases 
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provide  a  more  nearly  minimax  approximation  than  the 
one  obtained  through  use  of  the  methods  of  the  second 
class.  However,  most  algorithms  using  the  direct  approach 
may  require  many  iterations  to  arrive  at  a  good  approxi¬ 
mation,  especially  if  the  initial  approximation  is  poor. 

It  appears  that  an  investigation  of  the 
problem  of  which  combination  of  degrees  of  pn(x)  and 
q  (x)  -  where  n  and  m  are  variable,  but  their  sum  has 
an  upper  limit  -  leads  to  the  least  absolute  value  of 
minp  and  an  approximation  that  has  no  poles  in  the  interval 
of  approximation  would  be  of  theoretical  interest  and 
considerable  practical  value . 
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APPENDIX 


Remez’s  Second  Algorithm  adopted  to  Rational  Approximation 


0. 

1. 

COMMENT : 


2. 

3. 

4. 

5. 

6 . 

7- 

8. 

9. 

10. 

11. 

12. 

13. 

14. 

15. 

16. 
17. 


Steps  2  to  23  set  up  the  system  of  equations  for 
p^  and  . 

i  < —  1 

y(i)  -  (-1)1'1  pU) 

if- 1+1 

i  :  n+rrri-2., _ = _ >  3 

i  4 —  1 

j  < — 2 

M(i,  J-l)  ^-Y(i)[(x.)  J'"1] 
j  <—  j+1 

j  :  m+±,  — =— ; >  8 
i  < — i+1 

i  :  n+m+1,  _ ir _ ^  7 

i«—  1 
j  <—  m+2 

M(i,J-l)| - (x.)J''m'2 

j  <r~  j+1 

j  :  n+m+2,  — _ >  15 


18. 

19. 

20. 
21. 
22. 

23. 

24. 

25. 

26. 

27. 

28. 

29. 

30. 

31. 

32. 

33. 

34. 

35. 

36. 

COMMENT : 


A2  . 


1  <—1+1 

i  :  n+m+1,  —  ^  l4 

1  1 

M  ( i ,  n+m+2 )  < - Y(l) 

i  <—  1+1 

1  :  n+m+1,  _ 21 

Solve  system  of  linear  equations  to 
obtain  p0,  ...,pn  and  (%  =  l) 

i+— 1 

R(k,l) < — p( i-l) 
i  < —  1+1 

1  :  n+1,  — = — ^  26 
s(k,  1)  <— 1 
1<—  2 

i  :  m+1,  >  35 

s(k,i)  +—  q(l-l) 

1  < —  1+1 
Go  to  31 
k  « —  k+1 

k  :  3,  (<,  =  ,>)— >  (2,  37,  57) 

The  following  section  of  the  algorithm  gives 

the  equations  for  F[p(k)];  ie . ,  it  evaluates 

n  i  m  i 

.  Z  p.x  .  , 0  ,  and  . Znq.x  ,  IO  ,  and  f(x  ,  , 0 

1=0*1  n+m+2  }  1=1  i  n+m+2  J  v  n+m+2 

and  the  products. 
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37. 

38. 


39. 

40 . 

41. 

42. 

43. 

44. 

45. 

46. 

47. 

48. 

49. 

COMMENT : 


50. 

51. 

52. 


53. 

54. 


55. 

56. 


k  < —  1 
i  < — 2 

ts (k* 1) < —  0 

ts(k, l)  < —  ts(k, l)  +  s(k,l) C^m+n+1)1"1] 

i  :  m+lj  _ 40 

ts(k,2)*-(-l)n+m+2  [X  +  ts(k,l)  ] 


ts (k, 3) <—  o 

ts(k,3)  * —  ts(k,3)  +  H(l,l) [ (xn+m+2)1"1] 

i  :  n+1,  ^  »  46 

ts(k,4) f (xn+m+1)ts(k,2)  -  ts (kj 3) 

The  equations  for  F[ p (k) ]  are  given  by 
p (k) ts (k, 2)  +  ts (k, 4)  .  The  following  finds 


p ( 3)  by  the  secant  method. 


k  <—  k+1 

k  :  2,  ,  38 

p(l)ts(l,2)  +  ts(l,4) 
p(2)ts(2J,2)  +  ts(2,4) 

P ( 3) <—  p(2)  +  [ p ( 2 )  -  p(l)]E2/(E2-E1) 
k«—  3 


Go  to  2 


Jf . 


t: 

57. 

58. 

T: 

59. 

T: 

60. 

6l. 

62. 

63. 

64. 

65. 

66 . 

67. 

68. 

69. 

70. 

71. 

72. 

73. 

74. 

75. 


The  following  section  searches  for  Tp  using 
the  p^q^  associated  with  p(3). 

C1  cl/2 

c-^  :  .015*  — — 60 

This  gives  c-^  a  lower  bound  of  .015. 

C;L  <—  .015 

Steps  60  to  75  give  us  a  new  x-^  if  one  exists. 


h  <— 

2c1^x2”x 

v  <e- 

-p(3) 

xz  < 

“X1 

xz  < 

—  xz+h 

W 

-f(xz)  - 

Q(xz) 

w  : 

,r  < 

v*  ==— > 

76 

xz 

—  xz+h 

vw  f  (xz)  - 

Q(xz) 

vw 

:  v,  - 

-+72 

v  —  w 
w  vw 
G-o  to  66 

v  :  p  (3) »  — ■ g-— >  75 

Use  inverse  quadratic  interpolation  to 
find  new  x-^. 

Go  to  76 

Use  inverse  linear  interpolation  to  find 

new  xn . 

1 
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COMMENT:  Steps  76  to  97  give  a  new  x 


if  one  exists. 


i-1 

76.  i  < —  3 

77.  h^~ci(xi+i  "  xi-i) 


78. 

v  < — 

P ( 3)  (■ 

-1) 

i 

79. 

xz  <- 

1 — 1 

•H 

1 

+ 

h 

80. 

w« — 

f(xz) 

- 

Q(xz) 

81. 

1  w  1 

:  1  V  1 

> 

(<,  =  ,>)— >  (82,  85,  89) 

82. 

xz  <- 

1 — 1 

•H 

* 

- 

h 

83. 

w « — 

■f (xz) 

- 

Q(xz) 

84. 

1  w  1 

:  |  v  | 

> 

(<,  =  ,>)  — >(86,  85,  88) 

85.  xl-l< —  (xi_i  +  xz)/2 

86.  i  4 —  i+1 

87.  C-o  to  97. 

88.  h  4 - h 

89 .  xz  < — xz+h 

90.  vw  < — f(xz)  -  Q(xz) 

91.  |  vw  |  :  |  v  |  ,  _ £ _ *  95 

92 .  v  « —  w 

93  •  w  vw 

94.  Go  to  89 

95-  Use  quadratic  interpolation  (with  xz, 

xz-h3  xz-2h  and  vw,  w,  v  respectively  to 
find  new  x^  ^) . 

96 .  i  <r —  i+l 

97.  i  :  n+m+2,  _  —  77 


A6 . 


110 

111 

112 


00 

o\ 

h  «- 

-^cl(Xn+m+2 

xn+m+l ^ 

99. 

v  <- 

p(3)(-l)n+m+1 

100. 

xze-x  ,  .  n  +  h 

n+m+2 

101. 

W 

-f(xz)  -  Q(xz) 

102  . 

|  w 

V| 

> 

113 

103. 

xz  < 

r —  XZ+h 

104. 

vw< 

: - f (XZ)  -  Q(XZ) 

105. 

1  vw 

|  :  |  v  | 

109 

VO 

o 

1 — 1 

v  4- 

-  w 

107. 

w  <— 

-  vw 

108. 

Go 

to  103 

109. 

1  v 

1  :  p(3),  112 

Use  inverse  quadratic  interpolation  to 

find  new  x  .  , 0 

n+m+2 

Go  to  113 

Use  inverse  linear  interpolation  to  find 

new  x  ,  . 0 

n+m+2 


113. 

w  4 — 

f  (x1) 

1 - 1 

Gf 

1 

114. 

i  4— 

2 

115. 

V4 — 

f  (x±) 

-  Q(x±) 

116. 

1  V  1 

:  |  w  | 

< 

> - 

117. 

W  4— 

-v 

118. 

i  <■ — 

i+1 

119. 

i  : 

n+m+2 , 

< 

— — — > 

120. 

1  v  | 

:  P ( 3) 

+ 

123 


1 

:  rx 

A7. 


121.  Print  v,  p(i),  q(i) 

122.  STOP 

123.  p(l)  * — P ( 3 ) 

124.  (5/4)  p ( 3 )  :  I  V  I  ,  — 127 

125.  p(2) < —  |  v  | 

126.  Go  to  1 

127.  P (2)  (5/4) p ( 3) 

128.  Go  to  1 

129 .  End 


